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Nonlinear  Stability  of  Unsteady  Viscous  Flow,  Final  Technical  Report,  A.P.  Rothmayer 

The  stability  and  development  of  unsteady  separation  on  airfoil  leading  edges  has  been  investi¬ 
gated.  In  particular,  attention  was  focused  on  investigations  of  2D  unsteady  incompressible 
flow  past  a  parabola  at  angle  of  attack  which  models  the  leading  edge  of  many  airfoils.  This 
study  has  direct  application  to  leading  edge  stall  (LES)  and  thin  airfoil  stall  as  well  as  unsteady 
flow  past  leading  edges.  The  general  character  of  much  of  this  study  allows  for  application  to 
other  classes  of  unsteady  boundary  layer  flow.  The  accomplishments  are  detailed  below. 

1.0  Investigations  of  Boundary  Laver  Approaches: 

A  number  of  issues  were  examined  related  to  unsteady  leading  edge  separation  s<3lved  using 
the  classical  boundary  layer  approximation  at  high  Reynolds  numbers  (see  Fig.  1),  the  equa¬ 
tions  being: 


and 


Ux  +  Vy  =  0 


Ut  -I-  UUx  +  VUy  =  UCt  +  UeUCx  -f  Uyy  , 

where  Ue  is  the  inviscid  slip  velocity  at  the  solid  surface,  which  in  general  is  a  function  of  time 
and  position  along  the  solid  surface  and  is  specified  for  the  problem  in  question.  The  boundar 
conditions  are 


u(x,  0.  t)  =  v(x.O,  t)  =  0 


and 


u(x,  Y,  t)  -»  Ue(x,  t)  as  Y  -*  eo 
For  the  parabola  at  angle  of  attack  the  inviscid  solution  is 


Ue(x,t)  = 


y(x)  +  K(t) 

yy(x)2  -f  1 


and  y(s)  is  given  by 


yjy~  +  1  +  ln(y  +  /y-  -f  l) 
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The  angle  of  attack  parameter,  K(t),  measures  the  height  of  the  stagnation  point  below  the 
horizontal  axis,  and  may  be  related  to  the  physical  angle  of  attack  in  thin  airfoil  theory  when 
the  parabola  is  the  leading  edge  correction  for  flow  past  a  thin  airfoil.  The  K(t)  is  arbitrary. 
For  the  case  of  uniform  pitch— up,  K(t)  has  been  chosen  as 


K(.)  =  f 


,  1 ,  cosh(bt  -I-  c) 
b  cosh(c) 


1 


where,  a,  b  and  c  are  constants  and  a  is  the  ultimate  slope  of  the  K(t)  curve  (see  Fig.  2). 


a.)  The  existence  of  boundary  layer  instability  modes  for  flow  past  a  pitching  parabola  at  angle 
of  attack  was  verified  by  local  stability  analysis  (i.e.  Orr— Sommerfeld  analysis  applied  to  the 
boundary  layer  equations)  (see  Figs.  3  and  4).  Modeshapes  and  growth  rates  were  successfully 
compared  with  the  asymptotic  theory  of  Cowley,  Hocking  &  Tutty  (1985),  CHT,  (see  Figs.  4 
and  5).  The  stability  equations  were  developed  from  a  linear  perturbation  about  a  known  solu¬ 
tion  of  the  classical  boundary  layer  equations,  (uq,  Vq).  The  perturbation  takes  the  form 

(u,v)  ~  (uo,Vo)  -f  e(ui,vi)  +  ...  ,  e  <1  1  , 

1 

which  gives  the  perturbation  equations 
and 


UqUi^  UiUq^  +  VqUiy  +  =  ^lYY  ' 

These  equations  are  nonparallel  and  have  a  disturbance  growth  evolving  on  the  time  scale  of 
the  original  boundary  layer.  To  examine  high  frequency  instabilities,  we  made  a  partially  paral¬ 
lel  flow  approximation: 

(ui,vi)  =  e*«(u',v')  , 

which  gives  the  equations: 


iau'  -r  v'y  =  0 

u't  +  iauflu'x  -f  u'  +  Vqu'y  +  v'  Uqy  =  u'yy  • 
The  boundary  conditions  are: ' 


u'(x,0,  t)  =  v'(x,  0,  t)  =  0  . 


and 


u'(x,  co,t) ->  0  . 

This  system  was  solved  as  a  time  marching  problem  with  central  differences  in  Y.  The  growth 
rates  and  mode  shapes  were  computed  numerically  from  the  long-time  solutions.  The  above 
equations  are  a  non— asymptotic  high  frequency,  short  wavelength,  version  of  the  Cowley  et 
al  (1985)  work.  The  actual  asymptotic  structure  of  the  high  frequency  instability  is  a  bit  more 
complicated  than  this.  Any  high  frequency  instability  which  has  streamwise  wavelength 
will  be  inviscid  at  leading  order,  and  neutral  for  the  classical  boundary  layer.  The  CHT  instabil¬ 
ity  is  driven  by  viscous  effects  which  re-enter  through  a  viscous  critical  layer  centered  at  the 
minimum  velocity  point  in  an  unsteady  separation  or  maximum  in  a  jet  (see  Fig.  3). 

These  instabdities  should  render  the  unsteady  boundary  layer  equations  ill  -  posed  in  time  due 


to  the  fact  that  the  instabilities  at  zero  wavelength  have  unbounded  growth  rates.  However, 
the  instabilities  do  have  low  growth  rates,  and  so  a  smooth  enough  boundary  layer  solution 
which  evolves  more  rapidly  than  the  slow  growth  of  the  instabilities  should  be  computable. 

b. )  Joint  work  with  F.T.  Smith  concluded  that  the  above  CHT  instability  could  be  connected 
to  marginal  separation  in  a  limit  as  initial  separation  was  approached.  That  is,  the  unstable 
modes  with  infinite  growth  rate  at  zero  wavelength  continued  to  exist  until  the  point  of  first 
separation  (with  the  scaled  growth  rates  tending  to  zero  in  a  limit  as  the  first  separation  point 
was  approached). 

c. )  In  the  early  computations,  grid— grid  numerical  oscillations  were  encountered  whose  mo- 
deshapes  (see  Figs.  6  and  7),  growth  rates  and  overall  qualitative/quantitative  properties  coin¬ 
cided  quite  well  with  those  of  the  linear  boundary  layer  instability  computed  above. 


d.)  The  boundary  layer  oscillations  were  later  removed  in  the  full  parabola  at  angle  of  attack 
computations  by  an  appropriate  choice  of  numerical  scheme,  which  was  a  Crank -Nicholson 
method  with  the  farfield  conditions  removed  to  true  infinity  (see  Fig.  8).  We  believe  that  suffi¬ 
cient  smoothing  of  the  scheme  and  boundary  conditions  removed  the  oscillations. 


The  current  boundary  layer  equations  were  solved  in  streamfunction- velocity  form  with 
stretching  transformations  in  the  streamwise  and  normal  directions: 


s 


sinh(as) 

sinh(a) 


-  1  ^  S  ^  1  ,  -  Smax  ^  S  ^  S 


max 


N  =  N 


max 


sinh(bN) 

sinh(b) 


The  boundary  layer  equations  are  given  by 

u  =  0 


0  ^  N  ^  1  , 


0  ^  N  ^  N 


max  • 


and 


Ut  -f  SsUU^  -  SsNj^^Ui;^  =  Uct  -f  UeUcs  -b  N^nUj^^  -f-  NjliUjq^  . 


These  equations  were  solved  with  a  second  order  backward  temporal  difference  and  central 
differences  on  all  other  terms  including  the  streamwise  convective  terms.  The  method  is  glob¬ 
ally  iterated  to  convergence  at  each  time  level  with  each  alternate  sweep  being  in  opposite 
directions  to  accelerate  convergence. 


e.)  The  discrepancy  between  the  above  two  results,  i.e.  the  fact  that  an  Orr— Sommerfeld  lin¬ 
ear  stability  analysis  of  the  smooth  boundary  layer  solution  clearly  indicated  that  instability 
modes  should  be  present  and  would  be  dominated  by  the  shortest  wavelength  modes,  whereas 
the  numerical  computations  seemed  to  indicate  that  they  could  be  removed,  led  us  to  examine 
extensions  into  the  nonlinear  regime  of  the  linear  asymptotic  CHT  boundary  layer  instability 

\ 
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theory. 

The  nonlinear  critical  layer  is  structured  as  follows. 

The  nonlinear  structure  corresponding  to  the  linear  work  of  Cowley  era/  (1985)  is  found,  after 
some  trial  and  error,  to  occur  when  the  disturbance  within  the  critical  layer  (see  Fig.  3)  rise 
to  the  level  of  the  locally  small_garabolic  contribution  of  the  mean  flow.  In  the  critical  layer, 
it  is  nominally  expected  that: 

(y  —  V..')  ^ 

u  ~  Uo(y)  +  eU  +  ...  ~  Uo(yc)  +  (y  -  yc)Uo^(yc)  +  — 2 — ^Oyy(yc)  +  -  +  sU  +  ... 
For  simplicity  the  following  variables  were  defined:  ’ 

Uo,  =  U„(yc)  Ui,  =  Uo  &c)  Uo,  =  Uo_^(yc)  . 

The  critical  layer  is  placed  at  a  velocity  minimum  or  maximum  and  so  it  was  assumed  that 

Ui,  =  U„(yc)  =  0 

y 

To  preserve  the  connection  with  linear  theory,  the  critical  layer  thickness  was  taken  to  be  the 
Cowley  et  al  (1985)  value: 

y  =  y’c  ,  a  ^  1  . 


With  these  assumptions,  the  streamwise  velocity  expansion  becomes 

a-i/2Y2. 


u  -  Uo(y)  -  eU  -f 


U 


Oc 


•Uoc  +  ...  "T  eU  "T 


Therefore  the  perturbation  within  the  critical  layer  becomes  comparable  to  the  dominant  par¬ 
abolic  portion  of  the  mean  flow  which  drives  the  linear  boundary  layer  instability  when  the  per¬ 
turbation  rises  to  the  still  small  value  s  =  a  “  This  may  be  shown  to  generate  a  vertical  ve¬ 
locity  of  0(0^/^^ )  within  the  critical  layer  which  provides  the  displacement  driving  the  main 
linearized  boundary  layer  flow.  Pressure  displacement  interaction  is  negligible  provided  the 
wavelengths  of  the  disturbances  are  sufficiently  long. 


In  the  main  boundary  layer  flow  y  is  0(1)  and  the  streamwise  length  scale  is  taken  to  be  a  small 
specified  value  (see  Region  I  of  Fig.  3): 

ax  "ax 


Consistency  within  the  critical  layer  requires  that  the  disturbance  be  approximately  convected 
with  the  flow  velocity  (as  is  the  case  in  the  linear  CHT  mode)  and  so 


A. 

at 


-  ac 


1/2  a 

aT  • 


The  first  term  is  the  local  convection  with  the  flow  velocity,  where  Cq  is  the  wavespeed,  which 


\ 
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is  the  velocity  at  the  reversed  flow  minimum 


^0  “  ^Oc  ~  Uo(yc)  • 

The  second  term  is  determined  by  bringing  unsteady  effects  into  play  within  the  critical  layer. 
The  outer  flow  is  then  forced  by  the  displacement  effect  from  the  critical  layer  and  is  found 
to  have  the  expansions: 

u  ~  UQ(y)  +  T  a~^/^U2  +  a~-^/‘^U3  4-  ... 

and 

V  ~  4-  a^/^2  d-  4-  ...  . 

Three  orders  of  magnitude  are  required  to  reproduce  the  linear  results  presented  in  Cowley 
et  al  (1985).  In  addition,  the  critical  layer  has  algebraic  decay  in  the  matching  with  the  main 
boundary  layer  and  the  higher  order  terms_are  useful  for  numerically  imposing  that  matching 
boundary'  condition.  Substitution  into  the  unsteady  boundary  layer  equations, 


Ux  4-  Vy  =  0 


and 


Uj  +  UUx  -r  VUy  =  —  Px(X,  t)  4-  Uyy 

gives  the  continuity  equations 


Uix  +  Viy  =  0  i= 1,2,3 

and  the  momentum  equations 

(Uo(y)  -  Co)uix  +  Uo(y)viy  =  0 


(Uo(y)  -  Co)U2x  +  Uo(y)V2^  =  -  UjUix  -  ViU 


ly 


(Uo(y)  -  Co)’^3x  Uo(y)V3y  -  -  UiT  ~  ''l“2y  ''^2^1  y 

These  equations  were  integrated  to  give  the  expansions  in  the  main  boundary  layer 
u  ~  UoCy)  +  a-VdA(X,T)Uo(y)  4-  a-V2rB(X,T)Uo(y)  +  ^A^Uo'Cy) 


+  a 


-3/4; 


C(X,T)Uo(y)  +  ABUo(y)  +  ^A^Uo'Cy) 


+ 
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V  ~  -  a^/’Ax(Uo(y)  -  g  -  a'/2[Bx(U„(y)  -  Co)  +  AAxUo(y)] 

-  a‘HAx(U„(y)  -  c„)  +  At  +  (AB)xUi(y)  +  i(A3)xUo(y) 


where  A=A(X,T),  B=B(X,T),  and  C=C(X,T)  are  unknown  displacement  functions.  Notice 
that  the  smallest  term  entering-into  the  momentum  equation  in  the  first  three  orders  is 
0(0^/"* )  which  is  larger  than  any  0(1)  effect  coming  in  from  variations  in  the  mean  flow.  Also 
notice  that  the  leading  order  v-velocity  is  set  by  matching  with  the  critical  layer,  but  the  rest 
of  the  terms  are  forced  by  nonlinear  interactions  in  the  momentum  equation  (i.e.  the  power 
series  expansion  is  set  once  the  leading  order  v— velocity  is  known). 

The  scales  of  the  nonlinear  critical  layer  are  set  by  a  viscous  balance  and  by  assuming  that  the 
perturbations  to  the  main  boundary  layer  flow  enters  at  the  same  order  of  magnitude  as  the 
parabolic  mean  flow.  It  turns  out  that  long  scale  variations  from  the  main  boundary  layer  flow 
also  enter  this  balance  at  the  same  order  of  magnitude.  Therefore,  two  sets  of  scales  are  operat¬ 
ing  in  the  boundary  layer,  one  on  the  scale  of  the  local  perturbation  and  the  other  on“the  scale 
of  the  main  boundary  layer: 

dx  ^dX  BXq 


A  =  _  acn—  -f  +  — 

at  ^dX  dT  ato  • 

The  vertical  scale  of  the  critical  layer  is  determined  by  a  convective -viscous  balance  and  is 
found  to  be: 

.  ay  aY  ■  ayo 

The  perturbations  to  the  original  boundary  layer  flow  are  on  the  local  scales  (X,Y,T),  whereas 
the  original  boundary  layer  flow  is  a  known  function  of  the  longer  scales  (Xo,yo,to  )>  i-®- 
Uq  =  Uq  (xo,yo,  to)  •  The  expansions  in  the  critical  layer  are  found  to  be: 

U  -  Uo,(x„,yc.t„)  +  a-V2teu„,^^(x„.y„to)  +  U(X,Y,T)1  +  +  0-IU3  +  ... 


V  ~  aV4v(X,Y,T)  +  V2  +  a-V4v3  +  ...  . 

The  higher  order  terms  drop  from  the  equations,  and  are  set  from  matching  with  the  main 
boundary  layer.  Substitution  into  the  mass  conservation  equation  gives: 

Ux  4-  Vy  =  0 

The  momentum  equation  yields  the  following  equations  for  the  first  three  terms  in  the  expan- 


6 


Sion: 


(Uo(xo,yoto)  -  Co)Ux  =  0 
(Uo(xo,yc,to)  -  Co)U2x  =  0 

(UqCXq,  yc>  to)  t^o)^3x  ” ^Oy oy  1,(^0’ Vc’ to)  U  Ux  +  ^tUoy^^(xo,  yc,  to)V 

+  VUy  =  Uyy  -  P0x(X,  t)  +  Uoy^^  -  Uo  -  UqUoXo  ■ 

It  is  assumed  that  the  critical  layer  lies  at  a  local  minimum  or  maximum  in  the  velocity  profile 
and  so 

Uoc  =  Uoy/Xo,yoto)  =  0 

Furthermore,  the  disturbance  is  convected  with  the  local  flow  velocity  in  the  critical  layer 
which  implies  that 

^0  “  Uo(xo,yoto)  . 

At  the  critical  layer,  the  original  boundary  layer  satisfies  the  momentum  equation 

Uo,.  +  U„Uo^=  +  ■ 

Therefore,  the  first  two  equations  from  the  momentum  equation  are  identically  satisfied  and 
the  third  momentum  equation  becomes,  essentially,  a  classical  unsteady  boundary  layer  prob¬ 
lem: 

Uj  +  ^-Uoc  +  U  Ux  +  YU0(;V  -f  VUy  =  Uyy 

where 

Uoc  =  Uoy^o(%yc>to) 

is  a  known  constant.  Matching  the  critical  layer  with  the  main  boundary  layer  gives 

U(X,  Y,  T)  YAUo'c  +  ^Uo'c  as  Y  ^ 

and 

V(X,  Y,T)  ^  +  AY  +  ^A^  AxUoc  -  Aj  as  Y  «  . 

Within  the  critical  layer,  the  streamfunction  is 


7 


+  ... 


where 


il;  ~  a-^/'^Uo.Y  + 


Y3 


Uo'c  +  ^ 


U  =  V  =  -  » 


and  the  vorticity  is 


CO  =  Uy  —  Vx  ~  +  ... 


where 


^  “  ^xx  • 

From  the  above  discussion,  the  nonlinear  critical  layer  at  the  point  of  minimum  velocity  in  the 
reversed  flow  region  was  found  to  be  governed  by  the  following  system  of  equations,  obtained 
after  scaling  out  the  various  constants: 


Ux  +  Vy  =  0 

and 

U-j  ^^x  ^^FJy  ~  1  ^YY 

with  the  boundary’  conditions 


U(X,Y,T)-^+  YA  +  ^ 

as  Y 

V(X,  Y,T)  -  -  ^  +  AY  +  Ax  -  At 

as  Y  “ 

U(X,Y,T)->^ 

as  Y->  - 

V(X,Y,T)->0 

as  Y->  - 

Numerical  computations  of  the  early  initiation  phase  of  a  boundary  layer  separation  were  de¬ 
veloped  and  compared  with  existing  computations.  This  was  done  for  an  inviscid  vortex  travel¬ 
ing  over  a  plate  which  induces  an  unsteady  separation  in  the  boundary  layer  on  the  plate  (see 
Fig.  9).  These  computations  are  representative  of  current  state— of— the— art  boundary  layer 
computational  methods  in  the  early  pre— singularity  unsteady  separation  stage.  This  numeri¬ 
cal  method  is  also  the  same  method  used  to  obtain  smooth  solutions  on  the  parabola  at  angle 
of  attack  and  is  the  method  used  by  Peridier  et  al  (1991)  and  others  to  compute  the  early  stages 
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of  unsteady  laminar  boundary  layer  separation.  Wall  shears  for  the  vortex  convection  problem 
are  shown  in  Fig.  10  and  comparisons  with  Peridiereru/  (1991)  are  shown  in  Fig.  11.  Numerical 
computations  were  attempted  for  the  nonlinear  critical  layer,  described  above,  using  this  algo¬ 
rithm  (which  was  a  trivial  change  in  the  algorithm  from  the  above-mentioned  test  case).  The 
changes  were  due  to  the  coupling  algorithm  for  the  farfield  boundary  conditions.  The  finite 
difference  equations  are: 


1 

1 

X 

2 

yn+i  —  +  i  +  i 

_Jj _ y-t  ,  i-i,j  i-ij-i  _ 

2AY  2AY 


where  the  streamwise  derivatives  are  central  differenced  about  i— 1/2.  The  raoment'um  equa¬ 
tion  is  central  differenced  about  n + 1/2, j  after  substitution  of  the  mass  equation  for  the  stream- 
wise  convective  term: 


UP.+  i  -  U9. 

ij  ij 

’U'l+i  -t-  U"1 

+  Y^] 

AT 

2 

J 

2 

J 

'yn  +  l  +  vn] 

'U! 

2  J 

=  -  1  + 


Un  +  l  4-  U" 
^YY  YY 


The  resulting  equation  is  NeN^ton  linearized  and  coupled  to  the  mass  equation  using  a  block 
tri-diagonal  algorithm.  The  lower  boundary  conditions  are 


U"r‘  = 


and 


=  0  . 


The  upper  boundary  conditions  are  the  mass  equation  applied  at  j=N  and  the  central  differ¬ 
enced  form  of  the  matching  conditions.  The  v— matching  condition  is 


Vn  +  l/2 
i,N-l/2 


■yn  +  1/2  Y’  n  +  1/2 

i,N-l/2 


Y 


N-l/2 


^  pJi+l/2 
1 


AP  +  i  _  AP 
1  1 


AT 


This  equation  is  Newton  linearized  and  central  differenced  about  the  N— 1/2  gridpoint,  which 
is  taken  at  a  large  value  of  Y.  The  u— matching  condition  gives 


Un  +  l 


Y? 


(a-») 


where  g  denotes  guessed  values  from  the  Newton  linearization  which  are  iterated  to  conver¬ 
gence  at  each  time  level.  The  boundary  condition  is  implemented  in  the  quasi- simultaneous 
manner  of  Davis  &  Werle  (1985).  Following  Davis  &  Werle  (1985),  the  displacement  function 
computed  from  the  first  equation  (i.e.  A?^^ )  is  substituted  into  the  second  equation,  giving 
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a  single  boundary  condition  expressed  only  in  terms  of  variables  contained  within  the  critical 
layer  equations.  After  the  boundary  layer  equations  are  inverted  at  a  particular  i  location  the 
displacement  function  is  post— calculated  from: 


1 


Y2 

1  Tjn  +  l  N 

Y^  +  An+1  ‘’N  2 


21 


The  method  is  converged  at  each  time  level  by  repeatedly  sweeping  the  grid  in  X  from  a  fixed 
upstream  location  to  a  fixed  downstream  location.  The  method  is  converged  to  a  locql  absolute 
error  on  the  displacement  function  and  all  velocities  at  all  points  on  the  grid.  Stretching  is  used 
in  the  streamwise  direction  to  isolate  disturbances  near  the  origin  which  are  input  ^s  initial 
conditions  which  satisfy  mass  conservation. 

The  results  were  inconclusive.  Exponential  growth  was  observed  at  low  to  moderate  ampli¬ 
tude  (see  Figs.  12  through  15),  but  eventually  point-point  oscillations  occurred  with  no  clear 
way  to  remove  them.  Large  amplitude  limits  could  be  constructed,  but  without  a  numerical 
solution  to  connect  the  linear  instabilities  to  the  large  amplitude  limit  we  did  not  feel  that  we 
could  justify  the  limit  solution.  Over— clustering  the  grid  at  earlier  times  showed  that  these 
point —point  oscillations  could  be  moved  forward  in  time.  This  was  done  using  stretchings  com¬ 
parable  to  those  used  in  the  parabola  at  angle  of  attack  problem  and  the  vortex- over- a - 
plate  problem.  Our  feeling  at  this  point  is  that  there  is  a  problem  with  the  boundary  layer  equa¬ 
tions,  which  can  be  made  to  go  away  if  the  numerical  method  and  initial/boundary  conditions 
are  sufficiently  smooth  (and  the  above  critical  layer  is  under- resolved),  but  can  re-occur  if 
the  grid  spacing  is  made  small  enough. 

The  results  of  this  part  of  the  study  were  inconclusive  due  to  the  fact  that  we  cannot  obtain  a 
’’grid-independent”  verification  of  the  CHT  modes.  To  summarize 

1.  Non-  asymptotic  linear  stability  analysis  of  the  boundary  layer  equations  in  regions  of  flow 
separation  clearly  confirmed  the  existence  of  boundary  layer  instability  modes.  These  com¬ 
putations  are  smooth  and  regular. 

2.  Early  numerical  simulations  which  appeared  smooth  to  graphical  inspection  yielded  point - 
point  numerical  oscillations  which  were  in  complete  agreement  with  the  CHT  theory. 

3.  The  instability  encountered  in  #2  was  found  to  be  due  to  the  numerical  scheme  and  bound¬ 
ary  conditions.  It  is  believed  that  the  new  scheme,  which  now  conforms  to  current  state— of— 
the  art  boundary  layer  solution  strategies  for  early  separations,  is  smoother  than  our  first 
scheme.  The  boundary  layer  instability  modes  computed  from  linear  stability  analysis  in  #1 
are  still  present. 

4.  We  constructed  the  nonlinear  critical  layer  to  look  for  nonlinear  amplitude  modulation  and/ 
or  growth  continuation.  The  critical  layer  was  successfully  constructed,  but  the  numerical 
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method  used  in  #3  yielded  the  same  problem  encountered  in  #2.  It  was  found  that  the  numeri¬ 
cal  oscillations  were  controlled  by  grid  spacing  (which  was  the  same  conclusion  reached  in  #2). 
This  result  is  sensible.  Our  very  tentative  conclusion  is  that  boundary  layer  schemes  which  re¬ 
solve  the  critical  layer  will  eventually  run  into  numerical  difficulty  at  sufficiently  small  grid 
spacing.  This  suggestion  still  needs  to  be  verified. 


2.0  Stability  Analysis: 

Due  to  the  potential  difficulties  encountered  with  the  boundary  layer  equations  we  decided 
to  focus  part  of  our  attention  on  the  pre— separation  instabilities,  and  chose  the  Rayleigh  insta¬ 
bility  as  a  first  candidate  for  study.  Our  feeling  at  this  point  is  that  the  Rayleigh  instability  is 
one  of  a  number  of  possible  dominant  flow  solutions  which  may  occur  depending  on  the  fre¬ 
quency  of  the  laminar  flow  and  the  particular  flow  conditions  (The  likely  candidates  right  now 
for  2D  are:  Tollmien-Schlichting  waves,  laminar  boundary  layer  with  singularity  termination 
to  local  Euler  regions,  unsteady  marginal  separation  with  singularity  termination  to  local  Euler 
regions,  CHT  modes  -  which  are  likely  just  the  marginal  separations,  and  Rayleigh  instabili¬ 
ties).  We  believe  that  the  Rayleigh  instability  will  dominate  the  flow  near  a  leading  edge  sepa¬ 
ration  in  two  dimensions  providing  that  the  unsteady  forcing  (say  change  in  angle  of  attack) 
is  not  too  fast. 

a. )  We  performed  a  stability  analysis  for  the  primary  Rayleigh  modes  in  the  flow  past  a  pitching 
parabola  at  angle  of  attack,  and  verified  that  they  do  occur  prior  to  boundary  layer  separation. 
See  Fig.  16  for  the  neutral  stability  curve.  This  was  simply  a  solution  of  the  Rayleigh  equation 
given  a  boundar}'  layer  input  profile  and  is  substantially  the  same  as  the  other  Rayleigh  stability 
analyses  described  later  in  this  report. 

It  should  be  noted  that  the  eigensolutions  for  this  primary  boundary  layer  instability  do  pro¬ 
duce  eigenvalues  which  would  .allow  for  the  secondary  instability  cascading  discussed  below 
(see  Fig.  17). 

b.  It  was  recognized  that  a  series  of  cascading  secondary  instabilities  could  be  constructed  at 
low  disturbance  amplitude.  The  first  step  of  a  cascade  could  be  constructed  where  the  primary 
inviscid  instability  creates  a  Stokes  sublayer.  This  sublayer  becomes  unstable  to  secondary 
Rayleigh  instabOities  at  a  low  critical  disturbance  amplitude. 

The  classical  Rayleigh  instability  has  been  shown  by  Smith  &  Bodonyi  (1985)  to  occur  in  a  finite 
aspect  ratio  region  centered  within  a  classical  Prandtl  boundary  layer.  Long— wave  versions 
of  this  instability  for  triple- decks  have  been  considered  by  Tutty  &  Cowley  (1986).  The  scales 
for  the  dominant  Smith  &  Bodonyi  (1985)  instability  are  (see  Fig.  19,  Region  II): 

(x,y,t)  =  (xo,0,to)  +  Re-i/2(X,Y,T)  . 

Within  this  region,  the  streamwise  velocity  and  pressure  are  0(1)  to  match  with  the  oncoming 
boundary  layer  flow,  and  the  normal  velocity  is  finite  to  preserve  mass  conservation: 
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(u,v,p)~(U,V,P)  +  ... 


The  nonlinear  instability  is  governed  by  Euler  equations 

Ux  +  Vy  =  0  , 

Ut^  UUx  +  VUy  =  -  Px 

and 

Vt  4-  UVx  +  VVy  =  -  Py 

with  tangency  conditions  at  the  airfoil  surface  and  the  initial  conditions  being  the  rotational 
boundary  layer  flow.  In  streamfunction-vorticity  form  these  equations  become 

Q  =  T^xx  +  . 

and 

^  =  £2,  +  Wyax  - ’J'x£2y  =  0  . 

The  linear  version  of  this  problem  is  the  classical  Rayleigh  equation 

ia{Uo(Y)  -  c)[il)yy  -  a^]  -  iaUo"(Y)rp  =  0  , 

where 

Uo  =  'i^o'(Y)  , 

and  linear  normal —  mode  perturbations  have  been  assumed 

Q¥,Q)  ~  (^o(Y),Qq(Y))  +  s[(tp(Y),^(Y))ei“^^-='r)  +  c.c.'  . 

In  all  cases,  this  equation  was  solved  numerically  using  a  nonlinear  Newton  iteration  method, 
by  treating  the  complex  eigenvalue,  c,  as  an  additional  unknown,  with 

Cy  =  0  . 

Non -trivial  eigensolutions  are  enforced  by  requiring  that 

^y(0)  =  1  • 

Quasi— linearizing  and  central  differencing  of  the  above  equations  in  the  forms  indicated  al¬ 
low  inversion  of  the  resulting  finite— difference  system  as  a  set  of  coupled  block  tri— diagonal 
equations.  This  computation  becomes  difficult  near  neutral  points  due  to  the  logarithrmc  sin¬ 
gularity'  in  the  classical  Tollmien  expansions  about  the  critical  layer.  The  solutions  exhibited 
the  standard  properties  of  the  linear  Rayleigh  instability.  The  growth  rate  is  a  local  maximum 
within  the  Euler  region  and  is  stable  below  some  critical  wavelength.  The  disturbance  growth 
rates  also  become  small  as  the  scaled  wavelength  becomes  large,  which  connects  up  with  the 
long— wave  Rayleigh  instability  (see  Tutty  &  Cowley  (1986)  for  example).  For  dynamic  stall 
on  the  parabola  at  angle  of  attack  we  expect  all  inflectional  profiles  to  potentially  admit  Ray¬ 
leigh  instabilities  (though  this  is  not  a  sufficient  condition  for  instability).  The  unsteady  reverse 

V 
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flow  profiles  are  inflectional.  However,  boundary  layer  velocity  profiles  prior  to  flow  reversal 
are  also  inflectional.  Examination  of  the  numerical  boundary  layer  solutions  shows  that  the 
inflection  point  starts  at  the  wall  (as  opposed  to  the  flow  interior).  This  at  least  opens  up  the 
possibility  that  the  Rayleigh  instabilities  could  grow  out  of  a  Navier— Stokes  region  near  the 
wall.  However,  when  the  flow  near  the  wall  is  inflectional,  the  wall  shear  is  non— zero.  This 
means  that  the  velocity  profiles  in  any  near  wall  region  will  be  approximately  pure  shear  pro¬ 
files,  which  are  stable.  V^iatwefind  instead  is  that  the  neutral  solution  occurs  when  the  inflec¬ 
tion  point  is  at  a  finite  scaled  height  within  the  boundary  layer.  The  neutral  solution  is  an  invis- 
cid  neutral  mode  of  the  Smith  &  Bodonyi  (1985)  structure,  i.e.  the  classical  Rayleigh  problem, 
and  we  find  that  the  neutral  curve  for  the  dynamic  stall  occurs  somewhere  between  the  curve 
of  first  inflection  point  creation  (which  is  at  the  wall)  and  first  unsteady  flow  reversal. 

The  inviscid  Rayleigh  instability  does  not  satisfy  viscous  no  -  slip  conditions  at  the  wall.  There¬ 
fore,  a  viscous  sublayer  is  needed.  For  low  instability  amplitude  this  sublayer  will  be  a*  Stokes 
layer  (Region  III,  Fig.  19),  Can  this  Stokes  layer  be  destabilized?  The  answer  is  yes.  At  the  criti¬ 
cal  wave  amplitude 

u  ~  Uo(Y)  -f  s[u(Y)e‘“(^"‘='^^  -I-  c.c' 

with 


e  =  , 


a  Stokes  layer  is  generated  with 

(x,y,t)  =  (xo,0,to)  +  (Re-^/-X,Re-3/4Y,Re-i/2T) 

and 

(u, v, p)  ~  (Re - Re  “  Re “  ^/^P)  +  ...  , 
with  governing  equations 

iix  +  v.,.  =  0 

and 

Uj  =  —  Px(X,  T)  -f-  UyY  • 

Given  the  driving  disturbance  which  is  the  slip  velocity  of  the  linear  Rayleigh  instability  at  the 
edge  of  the  Stokes  sublayer,  the  exact  solution  is 

Ug  =  Uo'(0)Y  4-  u(0)(l  -  +  c.c  , 

where 

k+  = 

and  the  -b  indicates  that  the  complex  root  with  positive  real  part  is  chosen.The  wave  amplitude 
of  Re"^/"^  is  the  critical  case,  since  it  is  the  first  amplitude  at  which  a  Navier-Stokes  region 
can  be  created  within  the  Stokes  layer.  This  brings  nonlinear  terms  into  play,  which  will  allow 
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secondary  Rayleigh  instabilities  in  the  presence  of  sufficiently  inflectional  velocity  profiles. 
The  Navier- Stokes  region  has  the  scales  (see  Region  IV,  Fig.  19): 


(x,y,t)  =  (xo,0,to)  +  (Re-i/-Xo,0,0)  4-  (Re-^/^X.Re-V^Y.Re-iZ-T) 


and 

(u,v,p)  ~  (0,0, Re'^/^)  +  (Re-i/'^U,Re-i/4v,Re-^/2p)  +  _ 

The  region  is  governed  by  the  full  Navier— Stokes  equations: 

+  V  Y  =  0  , 

U’p  4"  UU^  4“  VUy  —  "1”  ^XX  ^YY 

and 

Vx  +  UVy  +  Wy  =  -  Py  4-  Vyx  +  Vyy  • 

We  solved  a  linear  version  of  the  above  equations  in  streamfunction  vorticity  form,  which  gives 
the  linear  stability,  Orr— Sommerfeld,  initial  value  problem: 

?  —  ’*i^YY  ” 

4-  iciUQ(Y,T)i5  —  icxij  qyy('^> ~  ^YY  ■ 

These  equations  cannot  be  solved  using  normal  modes  in  time,  due  to  the  fact  that  the  time 
scale  of  the  Stokes  layer  is  the  same  as  the  Navier-Stokes  region  and  T  appears  explicitly  in 
the  base  velocity.  This  linear  problem  was  solved  for  a  model  Rayleigh  instability  of  a  form 
similar  to  the  one  given  above,  but  with  arbitrarily  prescribed  growth  rates  and  other  constants. 
The  results  are  shown  in  Figs.  20  and  21,  where  the  average  streamwise  velocity  is 

Ly 

=  .  Ly  >  1  • 

0 

We  identified  two  typical  cases.  The  first  was  for  standing  wave  instabilities  (i.e.  zero  wa- 
vespeed  of  the  original  Rayleigh  instability)  in  which  the  solution  in  the  Stokes  layer  is  simply 
dragged  along  with  the  Rayleigh  instability  (see  Fig.  20).  The  second  case  was  a  traveling  wave 
instability  with  sufficiently  large  wavespeed  (see  Fig.  21).  In  this  case  the  solution  developed 
local  spikes  which  grew  at  a  faster  rate  than  the  main  Rayleigh  instability.  This  strongly  sug¬ 
gests  that  secondary  instabilities  could  occur  within  the  Stokes  layer,  and  could  develop  fast 
time  scales  and  large  amplitudes. 

b.)  We  developed  a  general  theory  of  cascading  linear  Rayleigh  instabilities  emerging  from  the 
above  structure  (see  Figs.  18  and  22).  We  showed  that  a  self— similar  discrete  scale  cascade  of 
instabilities  could  occur  which  asymptotes  from  the  boundary  layer  scale  to  the  viscous  dissipa¬ 
tion  scale  in  the  following  order  (for  linear  instabilities): 
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^n  + 1 


-  V  1/21^ 

=  Re 


k  =  l 


If  the  disturbance  amplitude  at  each  step  of  the  cascade  is  assumed  to  be  0(1)  and  it  is  assumed 
that  the  initial  instability  occurs  in  a  classical  boundary  layer,  then  the  cascade  scales  are: 


(Ai,A2,A3,A4,...)  =  (Re-i/2,Re-3/4^Re-7/8, Re-15/16,...) 


and  they  converge  quickly  to  the  dissipation  scale  Re  i .  Details  of  this  cascade  are  given  as 
follows. 

It  is  assumed  that  the  primary  Rayleigh  instability  occurs  in  a  region  of  dinlension 
(Aj  ,6^  =  A^ ),  where  for  the  classical  steady  or  unsteady  Prandtl  boundary  layer  (see  Fig.  23) 

=  Ai  =  Re-l/2  . 

Other  base  flows  may  be  used,  for  the:  natural  convection  boundary  layer  A-^  =  Re  “  i/'i ,  tri¬ 
ple-deck  separation  A  ^  =  Re-5/^,Tollmein-Schlicting  waves  Aj  =  Re “ ^/® ,  for  internal 
flows  separation  of  0(1)  streamwise  scale  A^  =  Re-^/^,  etc.  The  primary  instability  in  this 
first  step  of  the  cascade  is  governed  by  the  nonlinear  Euler  equations,  written  here  in  stream- 
function  vorticity  form 


^  =  Qt  -t-  UQv  +  ^^Y  =  0  • 

3)t  ^ 

For  all  the  Rayleigh  stability  calculations  in  this  study,  a  linear  stability  analysis  is  performed 
on  this  system  of  equations,  in  which  case 

CF,Q)  ~  (^o(Y).S2o(Y))  +  ei[(Tl)0/),?(Y))e'“^-"^)  +  c.c.] . 

Substitution  into  the  above  equations  gives 

i|jyy  -  ^ 


and 

ia(vFo'(Y)  -  c)|  -  iaWo'"(Y)rt»  =  0  , 

which  may  be  combined  to  give  the  classical  Rayleigh  equation 

ia(Uo(Y)  -  c)[rl)YY  “  a^]  “  iaUo"00^  =  0 
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where 


Uo  =  'i'o’m  ■ 

This  equation  may  be  solved  using  the  method  described  previously. 

The  inviscid  instability  of  0(8 1)  has  a  Xobem^^^^ 

Stokeslayerisneededtosatistyv~o-^^^^^^^^  ^  )_ 

perturbation  (i.e.  the  main  Rayleigh  instability)  are  given  by 

U  ~  Sn-l^n  +  - 

where  U„  is  the  base  velocity  which 


Note  that  for  the  classical  boundary  layer 


So  =  1 


..a «, ...  B S "• 

For  other  values  of  n  the  base  flow  is  a  .  inordinate  approaches  the  wall.  For  suffi- 

Rg-  24)  ,  .  V 


-  (x,y)  -  ■ 

We  assume,  for  now.  that  e,  - 

is  0(8n)  and  a  Stokes  layer  is  created  with  dimensions 


(x,y)~o^An,|^R^j  J  ■ 

A  new  secondary  Rayleigh  instability  can  occur  within  a  region  of  the  same  streamwise  leng 
as  the  height  of  the  Stokes  layer,  i.e.  ^  ^ 

(^’y)~^(R^)  ’(r^) 


and  so 


1/2 


•i+i  \  Res 


providing  that 


,Re  An 


En  ^  Sn_: 


The  lower  bound  on  e„  is  the  critical  perturbation  scale  for:  the  inviscid  repon  collapsing  into 
the  Stokes  layer  (in  which  case  the  Stokes  layer  is  driven  by  contributions  from  both  the  base 
flow  and  Rayleigh  instability),  and  that  the  secondary  region,  Region  V,  be  mviscid  If  En  is  at 
the  lower  bound  then  the  secondary  instability  region  is  viscous  and  governed  by  full  Navier 
Stokes  equations.  The  upper  bound  is  the  requirement  that  the  flow  m  the  P^mary  Rayleigh 
instability  be  linear  and  governed  by  the  classical  Rayleigh  equation  At  any  step  of  the  cascade 
it  may  be  easily  shown  that  the  secondary  Rayleigh  instability  problem  (i.e.  Region  V,  Fig.  24) 
is  given  by  the  classical  Rayleigh  equation 

ia(Uo(Y)  -  c)[^yy  -  “  iaUo”(Y)il)  =  0  , 

with  the  standard  boundary  conditions 

1|)(0)  =  \|)(co)  =  0  . 

The  base  flow  is  eiven  by  the  solution  to  the  boundary  layer  equations  in  the  first  step  and  by 
the  exact  solution  to  the  Stokes  layer  in  all  subsequent  steps,  which  is 


where 


Uo  =  u'n(Xo,0,To)(l  -  -b  c.c 


and  the  +  indicates  that  the  complex  root  with  positive  real  part  is  chosen,  see  Fig.  25. 

Various  properties  of  the  cascade  are  outlined  below.  As  noted  above,  bounds  on  allowable 
perturbations  at  each  step  of  the  cascade  may  be  found.  A  lower  bound  may  be  set  at 

"  (r^) 


and  it  is  assumed  that 


Efiov,  ^n-1 


It  may  be  easily  shown  that  if  En  > 
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1/2 


£n 

Sn-l 


6n  ^ 


An 


Res 


n-l 


which  means  that  the  inviscid  sublayer  (Region  III,  Fig.  24)  lies  above  the  Stokes  layer  (Region 
IV,  Fig.  24).  More  importantly,  if  Sn  ^  then  it  may  be  shown  that 


'■n+i 


<  tn 


which  means  that  the  time  scales  of  successive  Rayleigh  instabilities  is  decreasing  and  the  in¬ 
stability  develops  sooner  on  the  shorter  length  scales  which  are  closer  to  the  walk  Also,  if 
8n  >  Sn^  then  it  may  be  shown  that 

^n  +  l 


which  means  that  the  instability  cascade  occurs  on  successively  shorter  length  scales.  Finally, 
if  En  >  en„,  then  it  may  be  shown  that 


which  means  that  as  the  length  scale  of  the  secondary  instability  gets  smaller  it  takes  a  larger 
disturbance  amplitude  to  trigger  the  next  secondary  instability.The  recursive  expression  for 
the  scale  of  the  Rayleigh  instability  may  be  iterated  to  the  first  Rayleigh  instability  scale,  in 
which  case 


^n  +  l  “ 


-  y  1/2^ 
Re 


ii 

n 

k  =  l 


'n-k 


For  the  Prandtl  boundary  layer 

Eq  =  1  Ai  =  Re -1/2 


and  so 


^n  + 1 


-l/4n-  I  l/2'‘ 

Re 


n  —  1 

n 

k  =  l 


-1/2'^ 

n-k 


In  the  limit  as  n  =» ,  the  geometric  series  converges  to 


00 


and  so 
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^1  +  1 -"Re  ^  n  as  n  -  «  . 

k  =  l 

The  right  hand  side  is  bounded  by  the  inequality 

f]  =  Re-^'  • 

k=l  k=l 

Since  Sq  =  1  this  means  that  for  all  n 

An  >  Re~^  . 

« 

This  limit  is  achieved  rather  quickly  for  large  amplitude  disturbances.  For  example  if  the  dis¬ 
turbance  amplitude  at  each  step  of  the  cascade  is  assumed  to  be  0(1)  then  the  cascade  scales 
are 

(Aj,  A2,  A3,  A4, ...)  =  (Re"^/2,Re“^/^Re“'^/^Re~^^/^^,...) 

and  the  scales  asymptote  quickly  to  the  dissipation  scale  Re  “  ^ .  The  other  bound  is  easy  to  see 
from  the  individual  steps  of  the  cascade.  If  £n  =  £n^,  at  each  step  of  the  cascade  then  it  may 
be  easily  shown  that 


for  all  n.  The  scales  in  this  case  are  found  to  be 

A^  =  Re~^/^  ,  An  =  Re“-^/^  for  all  n  >  2  . 

The  physical  interpretation  of  this  result  is  simply  that  the  second  step  of  the  cascade  occurs 
in  a  foil  Navier— Stokes  region  and  so  no  more  cascading  is  possible.  This  does  form  a  reason¬ 
able  upper  bound  for  the  scales  of  the  cascade,  since  the  perturbation  at  each  level  may  be  tak¬ 
en  just  slightly  larger  than  the  critical  amplitude.  This  means  that  a  cascade  can  be  found 
where,  at  each  level,  the  spatial  scale  is  just  slightly  smaller  than  Re  Therefore,  the  cas¬ 
cade  scales  are  bounded  by 

Re"^  ^  An  Re-^/"^  . 

At  any  step  of  the  cascade  the  secondary  Rayleigh  instability  problem  (i.e.  Region  V,  Fig.  24) 
is  given  by  the  classical  Rayleigh  equation 

ian(Uo(Y)  -  CnfiprY  -  a^]  -  ia„Uo"(Y)4>  =  0 
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together  with  the  boundary  conditions 


a|,(0)  =  0  rp'(O)  =  1  =  0  • 

The  baseflow  is  simply  the  Stokes  layer  flow  generated  by  the  Rayleigh  instability  in  the  pre¬ 
vious  step  in  the  cascade,  which  is 

Uq  =  u'^(Xo,0,To)(l  -  4-  c.c 


This  equation  may  be  written  in  the  form 

Uo(Y)  =  Re‘®(l  -  e“^"^)  -f-  c.c  , 


where  0  is  a  parameter  that  positions  the  current  step  of  the  cascade  within  the  Stokes  layer 
and 


k+  =  ya„_i|c„_i| 


gi(Arg(c„_,)/2 -n:/4 +  (0,n)) 


and  the  +  indicates  that  the  complex  root  with  positive  real  part  is  chosen.  The  subscripts  n  - 1 
denote  values  which  are  known  from  the  Rayleigh  solution  in  the  previous  cascade  step  and 
the  subscript  n  denotes  values  in  the  current  step.  Note  that  the  eigenvalue  c  in  step  n  is  an 
unknown,  whereas  the  wavelength  a  in  step  n  is  a  specified  parameter.  Solving  the  above  prob¬ 
lem  directly  would  consist  of  solving  a  very'  large  number  of  Rayleigh  equations  coupled 
through  the  Stokes  layer  profiles.  At  each  step,  parametric  variations  in  streamwise  location 
within  the  Stokes  layer  as  well  as  wavelength  would  have  to  be  accounted  for.  It  is  clear  that 
only  2  or  3  steps  could  be  solved  by  a  direct  method.  Fortunately,  the  entire  cascade  problem 
can  be  reduced  to  an  iterative  mapping  on  a  reduced  Rayleigh  problem.  The  transformation 


(Y,tl;,a)  =  ^  {cn,Uo)  =  R(cn,Uo) 

V  *^n  — ll^n  — 1 1 

yields  the  following  reduced  problem;  the  Rayleigh  equation  and  boundary  conditions  are  the 
same  as  above,  but  the  base  velocity  is  replaced  by: 

Uo(Y)  =  e'®(l  -  -I- c.c  , 

where 

k+  = 

The  only  parameter  which  survives  from  the  initial  Rayleigh  instability  driving  the  Stokes  flow 
is  the  argument  of  the  complex  wavespeed  in  the  previous  step  of  the  cascade.  The  only  other 
two  independent  parameters  in  the  problem  are  ccn  and  0,  but  they  do  not  affect  the  next  step 
of  the  cascade.  Therefore,  the  only  thing  we  need  to  check  is  whether  or  not  the  flow  is  unstable 
and  then  compute  arg(cn)  in  the  unstable  regions.  An  interesting  consequence  of  this  trans¬ 
formation  is  that  the  entire  secondary  instability  cascade  is  completely  determined  by  the  argu¬ 
ment  of  the  eigenvalue  of  the  primary  instabOity  in  the  original  unsteady  boundary  layer  (i.e. 
by  arg(c)  computed  from  the  instability  of  the  original  boundary  layer).  A  schematic  diagram 
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of  the  reduced  problem  is  shown  in  Fig.  26. 


The  above  equations  are  Newton  linearized  and  finite  differenced,  which  gives 

Uo-cf 

+  a"|  Un  -  < 


Uo  -  cp 


AY2 


fj-i  - 


■^  +  a-=lU„-cf|  +  Uo 


+ 


Uo  -  cf] 


AY2 


^j  +  i 


+ 


C:  = 


cS( 


and 


[-  l]Cj_i  +  [l]Cj  =  0 


with  the  boundary  condition  y=0  at  the  wall, 

tpi  =  0 

and  the  second  order  accurate  boundary  condition  tpy  =  0  at  the  wall, 

=  AY 


The  free— stream  boundary  condition  is  imposed  by  considering  the  solution  to  Rayleigh’s 
equation  in  the  region  where  Uq  =  0 ,  this  gives  the  solution  il?  «  e  which  means  that  y 
satisfies  the  equation  'tpy  4-  aT|)  =  0.  This  last  equation  is  finite  differenced  in  a  second  order 
accurate  manner  to  give  the  boundary  condition; 


2. 

AY2 


o  1  +  aAY 
AY2 


a- 


=  0  • 


This  system  of  equations  is  inverted  using  a  complex  valued  block  tri-diagonal  algorithm.  A 
good  initial  guess  is  needed  for  this  algorithm  to  converge,  due  to  the  fact  that  there  are  nonuni¬ 
que  solutions  to  the  Rayleigh  equation.  In  general,  unstable  modes  coexist  with  stable  modes 
and  the  algorithm  has  a  tendency  to  lock  onto  the  stable  modes.  We  generate  the  initial  guess 
either  from  a  previously  converged  solution  or  by  solving  the  initial  value  problem  for  the  Ray¬ 
leigh  equation,  written  in  streamfunction  vorticity  form 

^YY  -  =  I 


and 


+  iaWo'(Y)^  -  iaT^o"'(Y)^  =  0  , 


with  boundary  conditions  that  mimic  a  fixed  wall  roughness: 
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tl)(0,T)  =  1 


These  equations  are  marched  long  enough  in  time  so  that  the  instability  modeshape  is  fully 
developed.  The  solution  of  this  initial  value  problem  is  only  needed  for  the  very  first  computa¬ 
tion.  Subsequent  eigenvalue  computations  use  initial  guesses  from  the  previously  computed 
eigensolution.  _ _ 

Since  the  only  parameter  from  the  n  - 1  step  driving  the  n’th  step  of  the  cascade  is  arg(c)  from 
the  n-1  step,  self  similar  solutions  occur  at  fixed  points  where  the  arg(c)  computed  at  the  n 
step  is  equal  to  the  arg(c)  input  as  a  parameter  from  the  n  - 1  step  (see  Fig.  26).  The  full  reduced 
problem  was  solved  near  a  region  where  self— similar  solutions  (i.e.  fixed  points)'were  ob¬ 
served  possible  (see  Figs.  27  and  28).  Self  similar  solutions  of  this  scale  cascade  were  computed 
and  are  shown  in  Figs.  29  and  30.  Fig.  27  shows  fixed  wavelength  slices  through  the  unstable 
region.  Fig.  29  shows  one  of  those  slices  with  contours  of  the  n— step  arg(c)  and  the  fixed  point 
surface  indicated.  Fig.  30  shows  the  function  arg(c)  along  the  cuts  1  -5  indicated  in  Fig.  29.  The 
passage  from  one  step  of  the  cascade  to  the  next  occurs  by  specifying  arg(c)  in  the  n— 1  step 
and  using  it  to  generate  arg(c)  in  the  n  step.  The  new  arg(c)  then  becomes  the  input  arg(c)  at 
n — 1  for  the  next  step  of  the  cascade  and  so  an  iterative  map  is  generated.  The  cascade  solutions 
shown  in  Fig.  30  assume  that  the  instabilities  in  the  various  steps  of  the  cascade  have  the  same 
scaled  wavelength  (i.e.  relative  to  the  Reynolds  number  scaling).  Furthermore,  it  was  found 
that  the  instabilities  with  maximum  growth  rate  at  each  step  of  the  cascade  could  also  form 
self— similar  cascades  to  the  dissipation  scale  (see  Fig.  28).  The  essential  idea  is  that  any  solu¬ 
tion  with  fixed  arg(c)  and  fixed  relative  position  in  the  Stokes  layer,  0,  will  have  a  single  maxi¬ 
mum  growth  rate  at  a  finite  value  of  wavenumber,  a.  This  means  that  the  waves  with  maximum 
growth  rate  form  a  surface  which  is  roughly  parallel  to  the  (arg(c),0)  plane.  If  that  surface  inter¬ 
sects  the  the  fixed  point  surface  in  such  a  way  that  contours  of  arg(c)  cut  across  the  fixed  point 
line  (i.e.  the  intersection  of  the  two  surfaces)  then  solutions  in  the  maximum  growth  rate  sur¬ 
face  look  similar  to  Fig.  29  and  cascading  solutions  can  be  found.  Fig.  28  shows  that  the  maxi¬ 
mum  growth  rate  surface,  mg,  does  indeed  intersect  the  fixed  point  surface,  fp.  Similar  scale 
cascades  are  likely  for  other  boundary  layer  instabilities.  The  one  which  dominates  is  likely  to 
depend  on  the  particular  solution  considered.  The  present  Rayleigh  instability  cascade  will 
likely  occur  for  loaded  airfoils  approaching  leading  edge  separation  and  lower  frequency  prob¬ 
lems  near  separation. 

It  is  clear  that  the  present  instability.  If  continued  to  the  viscous  dissipation  scale,  will  yield  a 
fractal  structure  in  the  streamwise  direction  at  infinite  Reynolds  number. 


3.0  Local  Numerical  Computations: 

The  local  nonlinear  problems  for  Rayleigh  instabilities  in  the  boundary  layer,  or  any  step  of 
the  cascade,  are  governed  by  Euler  equations  written  in  local  Cartesian  coordinates: 
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Ux  4-  Vy  =  0  , 

Uj  +  UUx  +  VUy  =  -  Px 


and 


Vt  +  UVx  +  VVy  =  -  Py 
In  streamfunction  vorticity  form  these  are 


^  ^YY  j 

^  =  Qt  +  UQx  +  VQy  =  0  . 

-  3)t 

A  model  problem  to  simulate  the  approach  of  the  above  structure  to  the  viscous  dissipation 
scale  is  the  Navier— Stokes  problem: 

Q  =  'PxX  -  'i^YY  . 

+  UQx  +  VQv  =  Re“i(Qxx  +  ^yv)  • 

A  wide  variety  of  numerical  methods  were  employed  to  compute  the  local  nonlinear  develop¬ 
ment  of  one  step  of  the  cascade.  For  the  sake  of  conciseness,  only  the  schemes  which  are  cur¬ 
rently  the  most  viable  will  be  discussed  below.  Convective  Taylor  series  expanded,  explicit  cen¬ 
tral  difference  and  implicit  central  difference  schemes  were  used  to  compute  both  Euler  and 
Navier -Stokes  equations  in  an  attempt  to  find  a  two  dimensional  soliton  emerging  from  the 
initial  linear  instability.  Local  Navier- Stokes  solutions  at  lower  Reynolds  numbers  did  show 
this  soliton  like  behavior  (see  Figs.  31  through  33) 

Three  stages  were  observed  in  this  growth  to  the  soliton  candidate.  The  first  was  the  classical 
cats— eye  pattern  with  linear  growth  about  the  critical  layer,  the  second  was  a  stage  with  a  main 
nonlinear  eddy  and  multiple  smaller  eddies  and  the  third  was  a  stage  with  a  single  large  scale 
eddy  of  finite  aspect  ratio.  The  creation  of  asymptotically  small  scale  eddies  is  believed  to  be 
unlikely  in  the  nonlinear  inviscid  equations  which  make  up  each  step  of  the  cascade,  due  to  the 
fact  that  all  functions  of  vortcity  must  be  conserved  in  an  inviscid  flow  which  tends  to  mitigate 
against  the  development  of  asymptotically  smaller  scale  structure  (note  that  many,  but  not  aU, 
cases  have  been  ruled  out). 

First,  we  note  that  the  full  nonlinear  coupling  in  the  Euler  equations  can  be  preserved  for  short 
wavelength  regions,  with  the  scalings 
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(X,Y,T)  =  6(X,Y,T) 


where  6  is  the  specified  small  spatial  extent  of  the  vortex.  The  dependent  variables  must  satisfy 

(T^.Q)  =  (6W,6“^Q)  . 


These  last  conditions  are  required  to  preserve  0(1)  velocities,  which  seems  reasonable  since 
the  wall  slip  velocity  and  all  other  velocities  within  most  of  the  small  scale  Euler  region  are 
likely  to  be  finite.  An  important  conclusion  is  that  the  vorticity  of  the  small  scale  vprtex  must 
become  locally  large  in  order  to  preserve  a  fully  rotational  flow.  Now  we  turn  to  the  conserva¬ 
tion  of  vorticity.  The  vorticity  equation  may  easily  be  manipulated  into  the  following  form 

=  (F(Q)),  +  T^y(F(^2))x  -  =  0  • 


In  vector  notation,  for  a  fixed  control  volume,  this  equation  is  equivalent  to: 


i. 

dt 


F(Q)dV 


F(Q)V  •  ndS  =  0 


where  V  is  a  fixed  volume,  S  is  the  bounding  surface  of  that  volume,  n  is  the  unit  outward  nor¬ 
mal  vector  to  the  surface  S  and  V  is  the  velocity  (Wy,  “  ^x)-  we  assume  that  the  volume  is 
a  periodic  box  with  one  surface  on  the  plate,  another  parallel  to  the  wall,  but  at  large  Y,  and 
the  period  taken  to  be  the  prescribed  period  of  the  flow,  then  it  may  be  easily  verified  that  the 
spatial  mean  of  any  function  of  the  vorticity  is  conserved.  Written  in  the  two-dimensional 
coordinate  system  this  means  that 


(F(Q))  = 


L/2  00 

I  F(Q)dYdX  =  const  . 

-L/2  0 


We  note  that  periodicity  is  not  likely  to  play  an  important  role  in  the  final  conclusions  drawn 
from  the  above  integral  equation.  For  example,  in  a  ’’localized”  disturbance  we  could  take  the 
edges  of  the  box  to  be  far  from  the  disturbance  and  arrive  at  the  same  result.  Such  a  formulation 
would  allow  for  complex  spatial  evolution,  but  would  be  difficult  to  verify  computationally. 
The  above  result  has  a  simple  physical  interpretation.  Each  differentially  small  patch  of  fluid 
has  a  well  defined  mean— vorticity  attached  to  it  (i.e.  approximately  equal  to  the  area  of  the 
fluid  patch  multiplied  by  the  magnitude  of  the  local  vorticity).  In  an  inviscid  flow,  the  vorticity 
is  convected  with  the  patch  of  fluid.  This  means  that  all  vorticity  within  the  flow  will  be  pre¬ 
served,  along  with  all  functions  that  can  be  formed  from  the  vorticity  (note  that  incompressibil¬ 
ity  and  two -dimensionality  is  an  important  part  of  this).  Particular  cases  that  were  examined 
in  the  numerical  solutions  were:  conservation  of  mean  vorticity,  mean  square  vorticity,  mean 
quartic  vorticity,  which  gives 
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L/2 


00 


(Q-)  = 


Q'^dYdX  =  const 


-L/2  0 


for  m=  1,2,  and  4,  as  well  as  mean  exponential  vorticity,  which  gives 

L/2  00 

(e^  -  l)dYdX  =  . 

Conservation  of  even  powers  of  vorticity  greater  than  2,  as  well  as  the  exponential*  vorticity, 
do  not  allow  for  the  creation  of  small  scale  structures  with  finite  velocity  that  are  governed  by 
the  full  Euler  equations.  This  is  because  of  the  fact  that  the  total  integral  over  suchj-egions 
would  yield  values  which  are  much  greater  than  those  present  in  the  original  flow  (which  is  set 
be  the  Stokes  sublayer).  It  is  also  believed  that  these  conservation  laws  will  play  a  key  role  in 
constructing  global  nonlinear  descriptions  of  the  dynamic  stall  boundary  layer  with  embedded 
eddies. 


-LjlO 


In  the  numerical  computations  it  is  important  to  maintain  both  stable,  oscillation  free,  com¬ 
putations  as  well  as  conservation  of  different  functions  of  vorticity  which  will  guarantee  that 
the  smaller  scale  structure  will  not  be  created. 


When  numerical  methods  used  for  the  Navier-Stokes  equations  were  applied  to  the  Euler 
equations  they  could  not  get  past  the  intermediate  nonlinear  multiple  -  eddy  stage  (note  that 
the  Euler  computations  are  at  very  high  Reynolds  numbers  and  are  the  computations  most 
relevant  to  the  scale  cascade  and  asymptotic  theory).  The  most  successful  straightforward 
computations  were  a  central  difference  method  in  conservation  form.  It  was  found  that  the 
central  difference  schemes  could  be  progressed  further  and  further  ahead  in  time  with  re¬ 
peated  spatial  grid  refinement.  However,  progress  was  slow  and  a  lot  of  spatial  grid  refinement 
was  needed  to  make  even  modest  gains  in  time.  The  finest  grids  run  in  our  computations  were 
500  by  500  spatial  grids  with  20,000  time  steps.  Central  difference  schemes  were  found  to  con¬ 
serve  vorticity  functions  well  (see  Fig.  38),  but  suffered  from  spurious  oscillations  within  the 
middle  of  the  multiple  eddy  stage.  Low  order  upwind  schemes  could  eliminate  the  numerical 
oscillations  and  pass  through  the  multiple  eddy  stage  to  the  later  stage,  but  failed  the  vorticity  - 
function  conservation  test  quite  badly  at  very  early  times. 

One  scheme  used  for  the  Euler  and  Navier— Stokes  solutions  consisted  of  convecting  the  vort¬ 
icity  and  then  Taylor  series  expanding  the  convected  vorticity  back  to  a  fixed  grid.  This  scheme 
was  used  for  both  the  inviscid  and  viscous  solutions.  It  was  found  that  a  conservative  central 
difference  scheme  yielded  substantially  the  same  results  as  the  following  scheme.  Both  the 
convective  and  central  difference  scheme  are  currently  being  abandoned  in  favor  of  the  CUD 
(compact  upwind  difference)  scheme  to  be  discussed  later.  The  three  schemes  are  outlined 
here  for  the  Euler  equations. 

For  the  Euler  equations,  a  velocity  field  is  calculated  at  the  n  - 1  time  step  by  central  differenc- 

\ 
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ing  the  streamfunction.  The  points  on  the  Eulerian  grid  are  convected  with  the  flow  via  an  ex¬ 
plicit  time  difference: 

Ax  =  +  O(AT-)  , 

-  W^-iAT  +  O(AT-)  . 

The  above  velocities  are  evaluated  with  second  order  accurate  central  differences.  The  above 
equations  can  be  extended  to  second  order  accuracte  expressions,  which  are: 

Ax  =  4- 

and 

Aj)  =  -  W^-iAT  -  +  +  0(AT3)  . 

The  vorticity  of  a  convected  gridpoint  is  taken  to  be  the  same  as  the  original  vorticity  at  the 
corresponding  Eulerian  gridpoint 


Q?,  = 

i  .  i 


The  vorticity  on  the  Lagrangian  grid  is  then  used  to  generate  new  values  of  vorticity  on  the  Eul¬ 
erian  grid  via  Taylor  series  expansion,  which  for  the  first  order  scheme  is: 


=  Q".  -f  AXijQx  +  Ay-n^,,  +  0{Ax,Ay-) 


These  equations  are  combined  into  the  following  Poisson— like  equation  for  the  vorticity 


i+i,j  i-i.j 


2AX 


Ay,j 


Ql'.  -  Q?.  , 

1.J  +  1  i,j  - 1 

2AY 


Qt 


where  x  is  a  pseudo— time  used  to  accelerate  and/or  stabilize  the  solution  within  each  physical 
time  step.  For  the  second  order  scheme  the  above  becomes 


Q?r‘  =  QP,  +  +  Ay.,-'’' 


10 


10 


2AK 


2AY 


+ 


A-2.Qn  .  _  2Qn.  +  Qn  Aw -  2Q"  -  Q" 

,  i+l,j  _ y _ i-l.j  ,  •'iQ  t.J  +  1 _ ^0 _ ^ 


AX2 


AY2 


-f  AxjjAyj  jQxY 


These  equations  are  solved  with  a  Peaceman— Rachford  ADI  method,  and  converged  in  the 
pseudo— time.  The  spatial  cross— derivative  in  the  second  order  method  is  iteratively  updated 
within  the  ADI  method.  After  the  vorticity  has  been  obtained,  the  central  difference  form  of 
the  streamfunction  equation  is  solved  using  an  ADI  method  with  specified  vorticity 
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mn  _  own  4.  Wn  Wn  _  +  W" 

i+l.j  ,  i-l.j  ip  i+l-j  _  on  ,  m 
AX2  AY2  ‘-i 

The  entire  numerical  method  is  formally  second  order  accurate  in  time  space.  Higher  order 
accurate  temporal  schemes  were  tested,  but  did  not  yield  significant  improvements  over  the 
above  scheme  on  fine  temporal-grids  (as  is  to  be  expected).  It  should  be  noted  that  this  method 
is  equivalent  to  a  central  difference  strongly  implicit  method  if  the  velocity  evaluations  are  tak¬ 
en  to  be  implicit  and  the  streamfunction  and  vorticity  equations  are  coupled. 

It  should  be  noted  that  the  above  scheme  was  used  to  generate  early  solutions  and  has  been 
abandoned  in  favor  of  central  difference  and  CUD  schemes. 

A  number  of  central  difference  schemes  were  tested.  Of  those  tested,  the  most  successful  ex¬ 
plicit  central  difference  scheme  was  found  to  be 


or 


r-.i  -  2'1’5  + 

Qn+  1  _  Qn-1 

(’f  =  0 


1  _  Qj.- 1  _  N”"),,.,  -  ^ 


2At 


2AX 


2AY 


Both  of  the  above  schemes,  as  well  as  the  current  CUD  scheme,  were  solved  with  periodicity 
conditions  of  the  form 


^?.j  =  ^NO 


and 


Q”  =  Q" 


These  periodicity  conditions  did  not  significantly  slow  the  numerical  method.  The  actual  peri¬ 
odicity  conditions  used  in  the  computations  were  applied  by  extending  the  grid  one  point  on 
either  side  and  iteratively  updating  the  extended  points  from  the  solution  on  the  interior  of  the 
grid.  The  actual  conditions  used  were 


mn  =:  \I/n 


wn  =  run 
^N  +  l.j  ^3.j 


and 

\ 

\ 
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where  1  and  N+ 1  are  the  added  points  and  N- 1  and  3  are  on  the  grid  interior.  These  periodici¬ 
ty  conditions  were  applied  to  all  local  schemes,  including  the  CUD— 3  scheme  discussed  below. 
Some  use  was  also  made  of  fixed  extrapolated  periodicity  conditions  of  the  form  (also  used  as 
initial  guesses  on  the  above)  _ _ 

nj  = 

«;.i  =  +  0(Ar) 

'*'no  =  2'!';-'  -  +  o(At2) 


and 


Q 


n 

N.j 


=  2Q 


n  —  1 

IJ 


However,  it  was  found  that  this  did  not  significantly  enhance  the  convergence  of  the  method. 
The  limiting  factor  in  overall  convergence  was,  and  still  is,  the  solution  of  the  Poisson  equation 
for  streamfunction. 


Discussion  was  initiated  with  the  United  Technologies  Corporation  to  interact  with  Professor 
Andrei  Tolstykh  of  the  Russian  academy  of  sciences  to  learn  the  CUD  methods  (Compact  Up¬ 
wind  Difference).  CUD  schemes  are  high  order  upwind  (i.e.  one-sided)  differencing  schemes 
that  maintain  3— point  computational  molecules.  The  third  order  CUD— 3  scheme  was  suc¬ 
cessfully  applied  to  the  Euler  equations  and  appears  capable  of  passing  through  to  the  soliton 
stage  while  maintaining  reasonably  good  conservation  of  vorticity  and  a  minimal  amount  of 
spurious  oscillations.  The  implementation  of  this  scheme  is  not  sufficiently  far  along  for  results 
to  be  presented.  However,  since  it  was  implemented  during  the  grant  period  and  since  it  will 
likely  form  the  core  of  future  work  the  basics  of  the  scheme  will  be  discussed  below.  The  meth¬ 
od  presented  below  is  first  order  accurate  in  time,  and  first  order  accurate  in  the  viscous  terms 
(to  be  discussed  later).  Second  order  temporal  and  viscous  schemes  are  straightforward  modi¬ 
fications  of  this  scheme,  but  do  not  retain  the  natural  simplicity  given  below. 


The  streamfunction  equation  is  solved  using  an  implicit  central  difference  ADI  with  pseudo 
time  stepping  within  the  iteration 


Ax2 


Ay- 


The  vorticity  transport  equation  was  solved  with  a  spatially  3rd  order  accurate  CUD  3 
scheme 


\ 
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Qt  +  (uQ)^  +  (vQ)y  =  0 


This  equation  is  particularly  appropriate  for  testing  the  CUD  schemes,  as  the  schemes  are  in¬ 
tended  to  use  upwind  differencing  at  high  order  to  model  convective  terms.  To  construct  the 
CUD -3  differencing,  consider  the  upwind  approximation  of  a  derivative 

—  f  =  ^ 

dx  • 

Normally  the  equation  is  evaluated  at  i  and  a  one  sided  difference  of  the  derivative  yields  first 
order  accuracy.  Third  order  accuracy  is  achieved  by  distributing  both  the  right  hand  side  and 
the  left  hand  side  of  the  above  equation  over  a  tri  -  diagonal  grid  using  the  difference  operators 

A.±fi  =  %^^+0(Ax3)  . 

where  the  standard  upwind  operators  (taken  in  x  and  y  directions)  are  ~ 

A,.,  =  o"+ij-()!;  ,  Ax_  =  ()”  -or-ij 


Ay+  -  Oij  5  Ay_  -  Oij  ()j_l 

and  the  CUD -3  operators  are 

Ax+  =  -  +  ^Oij  + 


Ay+  =  -  +  ^0"j  +  i 

Ay-  =  ^Oj-l  +  +  l  • 

These  operators  may  be  written  in  terms  of  a  switching  parameter  s.  For  example,  in  the  x— 
direction 


Ax(s) 


1  -f  s 
2 


0 


n 

i+lj 
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Ax(s) 


2  -  3s 


—  3s  a"  I  ^  I  2  +  3S/\n 

—  Oi-ij  +  3Uij  +  I2"^''i  +  ij 


J 


where  the  switching  parameter  gives 


-1,  backward  difference,  -  operator 
0,  central  difference 
+1,  forward  difference,  +  operator 


The  same  form  holds  for  the  y-direction.  The  forward  and  backward  differences  are  formally 
3rd  order  accurate,  while  the  central  difference  is  formally  4th  order  accurate.  To  enhance  sta¬ 
bility  it  is  sometimes  desirable  to  smoothly  switch  the  upwind  operators  between  a  forward  and 
backward  difference.  For  example,  this  avoids  the  situation  of  a  single  grid -point  oscillating 
from  forward  to  backward  difference.  This  was  accomplished  by  using  the  smooth  switching 
parameter,  given  in  the  x— direction  by 

s  =  —  tanh 


Suj 

1.35S  Au  +  lujl 


w'here  S  is  some  arbitrary  large  number  (S=20  was  used  in  our  computations)  and  Au  is  roughly 
the  velocity  half  range  about  u=0  where  the  smooth  switch  is  applied  (that  is,  when  u= Au  the 
switch  is  at  0.95).  As  mentioned  before,  if  s=  4- 1,-1  the  CUD -3  scheme  is  3rd  order  accurate 
and  upwinded.  If  s =0  the  CUD — 3  scheme  is  4th  order  accurate  and  central  differenced.  Other 
values  of  s  give  first  order  accuracy,  but  these  are  confined  to  thin  regions  where  the  velocity 
is  low. 


Consider  the  application  of  the  CUD-3  scheme  to  the  vorticity  transport  equation  in  two  di¬ 
mensions 


Qj  +  (uQ)x  +  (vQ)y  =  0  . 


The  first  order  temporal  CUD -3  scheme  was  taken  in  our  computations  to  be 
3Q"  -  4Q".-i  -b  Qn-2 


2At 


+  iAr'(s)A,(s)(Ln)”  +  iA,-‘(s)A,(s)(vQ)5  =  0  . 


In  approximately  factored  form,  this  scheme  becomes 


2At 


3Ax^^(^^  ^x(s)u; 


^  ^  ^Ay(s)v; 


-  -Q""- 
3  y 


This  scheme  is  formally  first  order  accurate  in  time,  but  may  easUy  be  extended  to  second  order 
accuracy  by  adding  a  grid  function  to  the  right  hand  side  of  the  equation.  The  factored  scheme 
is  solved  in  the  standard  fashion  by  splitting  it  into  two  steps 
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and 


’i  +  |^A,(s)-iA,(s)u;]  (Q);;  =  Iq;-!  - 


^  ^  ^Ay(s)v5 


(^)rj  =  Q? 


Each  step  of  the  CUD -3  algorithm  is  solved  by  multiplying  through  by  the  appropriate  opera¬ 
tor,  giving: 

A,(sXQ);  +  |^A,(s)(uQ);  =  AXs)(f  ns-‘  -  iQ?-2 


and 


A,(s)(Q)|;  +  ^A,(s)(vQ);  =  A,(s)Q5  . 

Now  consider  the  actual  implementation  of  each  step: 

I)  X— direction: 

Ax(s)(Q);  -  II^AxsKuQ);  =  Axs){fn;-' 


which  implies  that 


2  ~  3s  p>n 
12 


2  +  3s 
12 


Qr.ii  + 


2At 

s  -  1 

3Ax 

2 

(uQ)^ij  -  s(uQ); 


+  ^(uQ)^ij 


=  A,(s) 


Qn-l 

U 


—  Q"“- 


A  linearization  was  used  in  which  the  u  and  v  velocities  were  treated  as  known  guessed  values 
(in  our  problem  these  equations  were  iterated  with  an  ADI  solution  of  the  Poisson  equation 
for  streamfunction).  The  above  gives  the  tri— diagonal  system: 


2  -  3s 
12 


At(s-  1) 

3  Ax  ‘-d 


2sAt  n 
3Ax  y 


Q" 


-f 


'2  +  3s 
12 


4- 


At(s  +  1)  - 
3  Ax  ‘■'■d 


Ax(s) 


which  must  be  inverted  for  Q  at  all  i,j  gridpoints.  Note  that  the  right  hand  side  is  easily  com¬ 
puted  by  applying  the  A(s)  operator  to  the  known  vorticity  at  the  old  time  levels. 
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ID  y— direction: 


Ay(s)(Q)“j  +  ^A,(s)(vn);  =  Ay(s)Qn 


The  above  equations  imply  that 

2  —  3s on  I  2Qn  I  2  +  3s on 
12  ^  3^-ij  +  12 


+ 


2At 


+ 


IJ 

s  +  1 


(vQ)rj+i]  =  Ay(s)Qn 


Again,  assuming  that  the  v  velocities  are  treated  as  known  guessed  values,  the  above  gives  the 
tri- diagonal  system: 


2  -  3s 
12 


At(s  -  1) 
3  Ay 


Q3_i  + 


2  2sAt  n  O'* 

3  3Ay  y  y 


2  +  3s  ,  At(s  +  1) 

12  '  3Ay  y  +  i 


Ay(s)Q|; 


which  must  be  inverted  for  Q  at  all  i,j  gridpoints.  Note  that  the  right  hand  side  is  again  easily 
computed  by  applying  the  A(s)  operator  to  the  known  values  of  Q  computed  during  the  x-  in¬ 
version.  Also  note  that  the  switching  operator  s  is  different  for  the  x  and  y  inversions.  Care  must 
be  taken  to  consistently  apply  the  appropriate  switch  and  boundary  conditions  to  the  above 
equations. 


For  example  the  boundary  conditions  for  the  Euler  equations  are  given  by  flow  tangency,  or 

V  =  -  =  0 


The  vorticity  equation  evaluated  at  y=0  with  v— 0  is: 

-r  uQx  ~  0  • 

This  equation  may  be  written  in  a  conservation  form  which  mimics  the  interior  solution: 

Qt  +  (uQ)x  -  Qux  =  0  . 


The  CUD— 3  scheme  for  this  equation  is 


3Qn  -  4QP.-1  4-  Qn.-2 


+  iA,(s)-‘A,(s)(u£J)5  +  n?^[A,(s)  ‘A,(s)u;]  =  0  . 


2At  Ax ' 

This  equation  may  be  written  in  a  form  similar  to  the  interior  equations  by  multiplying  the 
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equation  by  2Dt/3  and  carrying  terms  to  be  treated  as  known  to  the  right  hand  side 


”3  +  |^Ax(5)-‘A,(s)(uQ);  =  |Q!-‘  -  -  n5|^[A,(s)-‘A,(s)u|’ 


Now,  a  grid  function  g  is  defined  such  that 


g  =  Ax(s)  ‘Ax(s)uJ 


The  grid  function  g  is  found  by  inverting  the  tri— diagonal  equation 


or,  in  expanded  form. 


where 


Ax(s)g  =  Ax(s)uJ 


2  3s  I  2  _n  I  2  ~1~  3^  —  a  — 

12  12  Si+ij  ’ 


A,(s)  =  -  so"-  + 


The  vorticity  transport  equation  at  the  wall  becomes 

”3  +  ^a,(s)-‘a.(s)(uQ)5  =  |n;-‘ 


hQ?r~ 


_  Qn  2At_r  n 
y  3Ax  L^yJ 


A.(s)£2;  +  |^A,(s)(uQ);  =  a,(s)[|q!-'  -  -  ^nrj8"j 


In  expanded  form  this  last  equation  is 


2  —  3s  t  2At  S  “  1  ,,n  on 
"l2“‘^3Ax  2 


2  _  2At.  nlQti 
[3  3Ax  yj  y 


2  +  3s 


,  2At  1  +  S  n  Ion  _  A  fc-'\[4nn-l  _  ion-2  _  2 At  ongii 

3M^^  y  3  y  3Ax  y^y 


Note  that  the  CUD -3  operator  on  the  right  hand  side  of  the  equation  is  simply 

A,(s)  =  2^o?-ij+fo;+2^or„j . 
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The  switching  parameter  may  be  the  simple  switch  from  backward  to  forward  difference  or  the 
more  complex  smooth  switch: 


s  =  —  tanh 


Su 


ij 


1.35SAu  +  U; 


4.0  Global  Numerical  Navier— Stokes  Computations: 


The  Navier- Stokes  equations  in  streamfunction  vorticity  form  were  solved  for  flow  past  a  pa¬ 
rabola  at  angle  of  attack  ^ 

A  series  of  full  Navier— Stokes  codes  have  been  written  for  flow  past  a  parabolic  leading  edge. 
The  main  code  used  was  a  pseudo- ADI  with  central  differencing.  This  algorithm  was  based 
on  a  globally  iterated  parabolized  Navier- Stokes  code  coupled  with  an  alternating  ADI  step 
for  the  streamfunction  and  vorticity  diffusion  (see  Davis  (1972)).  In  one  version  of  the  code 
the  dependent  variables  in  the  Navier- Stokes  equations  were  split  into  an  inviscid  and  a 
correction  siich  that  the  sum  of  the  two  yielded  the  correct  viscous  variables.  The  nonlinear 
’’perturbation”  equations  for  the  corrections  were  solved  by  the  above-mentioned  method. 
Details  of  the  various  codes  and  solutions  are  given  below. 


The  Navier-Stokes  equations  written  in  Cartesian  coordinates  and  non-dimensionalized 
with  the  leading  edge  radius  of  curvature  of  the  parabola  are  given  by 

-  Q  =  4-  Wyy 

and 

Qt  +  WyQx  -  WxQy  =  Re"^(Qxx  +  ^yy)  • 


These  equations  were  transformed  to  parabolic  coordinates  using 


X  = 


y  =  Iq 


which  gives  the  equations 

-  (|2  +  ,,2)£2  =  +  W,, 

and 

(|2  +  r]-)Q,  +  . 

The  boundary  conditions  are  no— slip  at  the  parabola  surface 

'Tf^.l.t)  =  'Tr,(^,l,t)  =  0 
and  an  asymptote  to  the  free -stream 

mw.t)  = 
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and 


The  inviscid  flow  is  given  by 


=  (I  +  K(t))(ii  -  1) 

and 


^inv(i>n.t)  =  0  , 

where  K(t)  is  the  angle  of  attack  parameter  discussed  in  the  boundary  layer  solution!  The  fol¬ 
lowing  grid  stretchings  were  added  to  cluster  points  in  the  boundary  layer  and  in  the  stream- 
wise  directions  near  regions  of  interest 


i  ^  =  11(11)  • 

Using  these  transformations,  the  Navier- Stokes  equations  become 

-  (r- +  ii^)Q  = 


and 


+  -n-jQt  -f 


The  stretching  in  the  normal  direction  to  the  parabola  is  used  to  resolve  the  boundary  layer 
and  sublayers  and  is  given  by 


11  -  1 


sinh(anil) 


11  max  -  1  sinh(an) 


0  ^  Ti  ^  1 
1  =  11  =  Tlr 


To  place  points  near  the  leading  edge,  a  simple  clustering  transformation  was  used 

^  _  sinh(as^)  -  1  ^  I  ^  1 

imax  sinh(as)  —  ^^ax  —  ^  —  ^max 


The  ”a”  parameters  give  the  degree  of  stretching  in  each  case.  To  cluster  points  in  a  specific 
region  of  interest  in  the  streamwise  direction,  the  grid  was  split  into  three  regions: 


Region  1:  -  Imax  =  ^  = 

Region  2;  ^  |  ^  §2 

Region  3:  ^2  -  ^ 


-  1  ^  ^ 

^  I  ^  I2 
I2  ^  ^  1 


The  stretching  transformation  in  Region  1  is: 
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sinh  a,  1  -  = 


I  4-1 


+  Ir 


i  +  i, 


max 


In  Region  2  a  uniform  grid  was  used: 


sinh(ai) 


I  =  -  ^i)  ■ 

S2  ~  ?1 


In  Region  3: 


I  -  I2 

- — Y"  =  sinh 

^max  S2 


1-g. 
}  ~  ^2> 


\ 


The  a’s  are  chosen  so  that  the  strearawise  metric  is  smooth  throughout  the  domain.  A  typical 
grid  is  shown  in  Fig.  40. 


The  finite  difference  scheme  is  a  central  difference  pseudo- ADI  based  on  global  iteration  of 
a  parabolized  vorticity  equation.  The  first  step  of  the  ADI  is 


W,  -  (|-  +  Ti2)Q"-1/2  = 


-1/2 


■1/: 


where  the  first  term  is  a  fictitious  temporal  term  added  to  stabilize  the  method  and  is  iterated 
to  zero  at  each  real  time  level.  The  vorticity  transport  equation  is 


All  solid  underlined  terms  (at  the  n+1/2  level)  were  treated  implicitly.  All  dashed  underlined 
terms  were  assumed  known  from  the  n  time  level.  All  terms  were  spatially  central  differenced. 
The  n+1/2  time  level  was  iterated  to  convergence  using  a  global  iteration  procedure  in  the 
streamwise  direction.  That  is,  the  equations  were  central  differenced,  Newton — linearized  and 
inverted  in  the  normal  direction  using  a  block  tri- diagonal  algorithm.  Repeated  sweeps  in  the 
streamwise  direction  were  used  to  converge  the  n+ 1/2  time  level.  The  direction  of  the  sweeps 
were  alternated  with  each  pass  to  enhance  information  propagation  over  the  upper  and  lower 
surfaces  of  the  parabola. 

The  second  step  of  the  ADI  method  consists  of  a  streamwise  inversion  of  the  vorticity  equation 
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(r-  +  r,’-)- 


Qn  +  l  _  Q"  +  1/2 


At/2 


+ 


li;n  +  l/2Qn  +  l/2  _  n;n  +  l/2Qn  +  l/2 

’T  I  in 


Again,  solid  underlined  terms  at  the  n+1  level  are  unknown  and  dashed  terms  at  the  n+1/2 
level  are  known.  This  means  that  the  above  equation  could  be  inverted  as  it  stands  with  a  single 
equation  tri- diagonal  algorithm  applied  to  each  normal  grid  line.  After  this  is  accomplished, 
the  streamfunction  is  found  from 


-  (r  + 


Again,  solid  underlined  terms  are  unknown  and  the  fictitious  temporal  term  is  iterated  to  zero 
at  the  n+1  time  level.  The  above  scheme  is  formally  second  order  accurate  in  space  and  time. 


The  boundary  conditions  at  downstream  infinity  on  the  upper  and  lower  surfaces  of  the  parab¬ 
ola  were  assumed  to  be  asymptotes  to  the  following  Blasius  solution 


^  =  ?F(ri) 


These  equations  were  used  to  develop  expressions  for  the  various  streamwise  derivatives  ap¬ 
pearing  in  the  above  equations  at  the  first  and  last  streamwise  stations,  which  are 


and 


= 


Q 


=  0 


2Q 


The  above  Navier- Stokes  code  has  been  verified  for  the  symmetric  steady  case  against  older 
well -documented  steady  symmetric  solutions  (see  Fig.  41) 


We  have  used  this  code  to  compute  a  number  cases: 

a.)  This  Navier -Stokes  code  was  used  to  compute  an  impulsively  started  parabola  at  angle  of 
attack.  An  initial  inviscid  flow  is  used  to  compute  a  viscous  flow  by  impulsively  applying  the 
no-slip  condition  at  the  surface  of  the  parabola.  This  reproduces  results  of  Reisenthel  (1994) 
(see  Figs.  42  and  43).  Secondary  separation  was  observed  (and  will  be  described  below), 
though  we  do  not  beUeve  that  this  is  tied  to  the  instability  cascading  discussed  previously  in  this 
report.  The  impulsive  start  case  is  likely  to  be  related  to  the  laminar  boundary  layer  singulari- 
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ties.  It  should  be  noted  that  we  could  not  obtain  grid  independence  in  this  particular  problem. 
We  believe  that  this  is  due  to  a  failure  to  represent  the  unsteady  Rayleigh  layer  at  the  impulsive 
start-up  (i.e.  the  grid  is  taken  to  be  fixed).  The  reason  this  was  done  was  to  keep  these  particu¬ 
lar  computations  as  close  as  possible  to  others  currently  in  the  literature,  which  also  use  a  poor 
initial  start-up. 

b.)  Using  this  Navier- Stokes  c^e,  we  have  shown  that  a  similar  structure  occurs  in  the  case 
of  a  rapid  unsteady  pitch-up  (see  Figs.  44  through  48).  The  same  multiple-eddy  structure 
occurs  in  the  boundary  layer  (see  Figs.  44  and  45).  Here,  we  were  able  to  come  close  to  grid 
independence  (see  Figs.  46,  47  and  48).  Furthermore,  multiple  embedded  eddies  were  ob¬ 
served  (see  Fig.  48) 

The  Navier -Stokes  code  used  here  is  slightly  modified  from  the  one  discussed  above.  It  was 
found  that  the  free -stream  boundary  condition  could  be  more  effectively  imposed  on  the  inte¬ 
rior  flow  if  the  solution  was  split  into  the  inviscid  part  plus  a  nonlinear  perturbation  which  gives 
the  true  viscous  solution.  That  is 

=  ri;  -K 


and 


Q  =  0)  +  =  to 

since 


The  original  Navier-Stokes  equations  were 

-  (|2  +  ,i2)£2  = 

and 

(^2  +  =  Re-^(Q.t  +  Qr,n)  • 

Substituting  the  above  equations  for  the  transformations  gives  the  perturbation  streamfunc- 
tion  equation 

due  to  the  fact  that  the  inviscid  solution  satisfies  the  Laplace  equation,  and  the  perturbation 
vorticity  transport  equation 

(^2  +  ^  =  Re  ^(c0||  -I-  tOT,Ti)  . 

These  last  two  equations  for  the  perturbations  were  solved  in  the  same  manner  as  the  full  sys¬ 
tem.  The  boundary  conditions  on  these  equations  become 
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=  0 


and 


at  the  wall,  and  the  asymptote'to’the  free -stream 

^(^>^max?t)  —  W(^,  Tlrnax’ 0  0  . 


c. )  We  have  also  shown  that  impulsively  changing  the  angle  of  attach  in  small  increments 
(which  will  lead  to  a  sequence  of  quasi-steady  solutions  in  stable  cases)  will  give  rapid  multi¬ 
ple  eddy  breakup  in  a  flow  which  is  immediately  post  separation.  We  think  that  this  may  be 
connected  either  to  marginal  separations. 

d. )  Applying  the  quasi -steady  impulsive  change  in  angle  of  attack  in  pre"- separation  cases 
produces  unstable  long-lived  unsteady  flow  with  multiple  eddies  (see  Fig.  49).  The  times  in 
these  figures  are  75  time  units  which  is  somewhat  longer  than  the  30  time  units  in  the  rapid 
pitch  -  up  case.  A  blow-up  of  the  local  wall  shear  near  the  eddies  is  shown  in  Fig.  50  and.vortic- 
ity’  contours  are  shown  in  Fig.  51.  It  should  be  noted  that  very  fine  grids  were  needed  to  obtain 
these  solutions.  A  linearized  boundary  condition  was  added  to  the  parabola  surface  to  enable 
us  to  input  surface  distortions  of  small  amplitude  but  arbitrary  functional  form.  The  eddy  struc¬ 
ture  is  strongly  influenced  by  small  changes  in  surface  geometry  perturbations  (see  Fig.  52). 
Vorticity  and  streamline  patterns  in  the  vicinity  of  the  eddies  are  shown  in  Figs.  53  and  54.  Note 
the  resemblance  of  the  eddy  in  Fig.  54  to  the  earlier  computed  eddies  in  the  local  solution  (see 
Figs.  31  and  34).  We  believe  that  these  last  computations  (i.e.  impulsive  change  in  angle  of  at¬ 
tack  below  separation)  are  candidates  for  Rayleigh  instabilities.  It  is  clear  that  we  can  rule  out 
all  separation  eddy  creation  mechanisms,  such  as  marginal  separations  and  boundary  layer  sin¬ 
gularities,  due  to  the  fact  that  the  flow  remains  attached  while  the  initial  instability  developes. 
It  is  not  clear  at  this  point  what  effect  of  the  impulsive  change  in  angle  of  attack  has  on  the  com¬ 
putations.  This  should  be  changed  to  a  smooth  variation  in  future  computations. 

Lastly,  a  CUD-3  Navier- Stokes  algorithm  based  on  the  pseudo -ADI  parabolized  vortcity 
iteration  was  written  and  verified  for  the  steady  symmetric  solution  (Fig.  55)  and  early  stages 
of  the  unsteady  pitch -up  (Fig.  56). 

The  unsteady  CUD-3  Navier-Stokes  algorithm  uses  the  stretched  equations  in  parabolic 
coordinates  written  in  conservation  form 

XP,  (^2  +  ^2)Q 

J  J 


and 


V 
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where 


a  fir-  + 

1  +  (WjrQt-  (w^q)  =  Re~^ 

r 

+ 

at  1  j  1 

\  / 

< 

^  J 

The  CUD  -  3  algorithm  is  applied  to  the  above  equations  in  a  form  similar  to  the  pseudo  -  ADI 
algorithm  discussed  above.  Nominally,  the  CUD— 3  algorithm  is  given  by: 


where  the  A’s  and  A’s  are  the  normal  CUD -3  operators  and  upwind  operators  respectively 
and  the  A’s  on  the  right  hand  side  are  central  difference  operators,  given  by 


The  algorithm  is  implemented  in  a  manner  similar  to  the  pseudo— ADI  algorithm  discussed 
earlier.  Since  this  work  has  only  progressed  through  the  very  preliminary  stages,  further  details 
will  be  omitted  here.  The  above  scheme  is  second  order  accurate  in  time,  third  order  accurate 
in  space  in  the  convective  terms  and  first  order  accurate  in  space  in  the  viscous  terms  (which 
appears  to  be  a  standard  accepted  limitation  of  past  CUD  work).  We  believe  that  the  CUD 
schemes  will  be  a  major  emphasis  of  future  work  and  attention  will  be  paid  to  increasing  the 
accuracy  of  the  viscous  terms.  It  is  clear  that  the  accuracy  of  both  the  unsteady  terms  in  the  fully 
implicit  method  as  well  as  the  viscous  terms  may  be  made  second  order  in  space  and  time  using 
methods  that  should  work  well  within  a  Newton  iteration.  Despite  the  nominally  low  accuracy, 
full  viscous  solutions  we  computed  showed  good  agreement  with  the  central  difference  scheme 
for  steady  and  early  unsteady  flows. 
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Fig  2.  Temporal  variation  of  K(t)  for  smooth  constant  rate 
pitch-up.  All  solutions  are  started  from  the  steady  solution 
for  a  flow  past  a  parabola  at  zero  angle  of  attack. 
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viscous  critical  layer  — 


Fig.  3.  The  high  frequency  boundary  layer  instability  structure  of 
Cowley,  Hocking  &  Tutty  (1985),  showing  the  viscous  critical 
layer  positioned  at  the  velocity  minimum. 
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Fig  5.  Comparison  of  Boundary  Layer  Stability  Analysis 
Mode  Shapes  with  CHT  Theory. 
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Fig.  7.  Comparisons  of  grid-grid  oscillation  mode 
shapes  with  boundary  layer  stability  mode  shapes. 
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Fig.  8.  Wall  shear  for  new  boundary  layer  computa¬ 
tions  without  oscillations. 
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Fig.  9.  Geometry  for  a  vortex  moving  over  a 
flat  plate. 


Fig.  10.  Scaled  wall  shear 
separation  computed  at  v 
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early  in  the  unsteady 
times. 
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Fig.  11.  Displacement  thickness  computed  at  various 
times,  the  time  interval  between  plots  is  0.1.  Also  shown 
are  comparisons  with  Peridier  et  al  (1991).  Results  were 
electronically  scanned  from  original  reprints  of  articles. 
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Fig.  12.  Maximum  absolute  displacement  function  and 
displacement  function  slope  showing  the  low  amplitude 
exponential  growth. 
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Fig.  13.  Displacement  function  at  various  times, 
showing  the  growth  of  the  displacement  function 
with  time. 
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Fig.  14.  Deviation  of  streamwise  velocity  from  parabolic 
profile  in  the  critical  layer  at  various  times,  at  a  typical 
streamwise  location. 
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Fig.  16.  The  Rayleigh  instability  neutral  curve  for  a  pa¬ 
rabola  in  smooth  pitch -up.  Also  shown  are  the  separa¬ 
tion  curve  and  the  curve  of  first  inflection  point  creation 
at  the  wall. 
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Fig.  17.  The  arg(c)  for  the  primary  Rayleigh  instability 
for  a  parabola  in  smooth  pitch— up. 


nonlinear  evolution  linear  instability 


Fig.  18.  Anticipated  structure  of  the  Rayleigh 
instability  cascade. 
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Fig.  19.  The  initial  2-step  cascade:  I-boundary 
layer;  Il-linear  Euler;  Ill-Stokes  layer;  IV  linear  or 
nonlinear  Navier-Stokes  region. 
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Fig.  20.  Response  of  the  Navier-Stokes  region  of  the  two- 
step  cascade  to  an  inviscid  standing  wave  instability. 
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Fig.  21.  Typical  secondary  instability  in  the  Navier- 
Stokes  region  of  the  two-step  cascade  generated  by 
an  inviscid  traveling-wave  instability. 
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Fig.  22.  Schematic  diagram  of  the  linear  instability  cascade. 
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Fig.  23.  The  initial  large  amplitude  2-step  cascade  in  a 
classical  external  boundary  layer:  I-classical  boundary 
layer;  Il-linear  Euler  region;  III-  passive  thin  layer  Euler; 
IV-Stokes  layer;  V  linear  or  nonlinear  Euler  region;  VI- 
Stokes  layer. 
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Fig.  24.  A  general  step  of  the  cascade:  I-  the  initial  Stokes 
layer;  Il-linear  Euler;  HI-  passive  thin  layer  Euler;  IV- 
secondary  Stokes  layer;  V  linear  or  nonlinear  Euler  region. 
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Fig.  26.  Schematic  diagram  of  the  reduced  stability 
problem  for  a  typical  cell  in  the  instability  cascade. 
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Fig.  28.  The  right  fixed  point  surface  showing  con¬ 
tours  of  arg(c),  as  well  as  the  maximum  growth 
rate  surface,  showing  the  intersection  of  the  maxi¬ 
mum  growth  rate  surface  with  the  fixed  point  sur¬ 
face.  \ 
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Fig.  29.  A  fixed  wavelength  cut  through  the 
unstable  region. 
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Fig.  30.  Arg(c)  along  the  cuts  of  Fig.  29  showing  an 
unstable  fixed  point  and  a  stable  fixed  point.  The 
stable  fixed  point  yields  a  self- similar  cascade  to 
the  viscous  dissipation  scale. 
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Fig.  31.  Streamfunction  contours  for  the  nonlinear  de¬ 
velopment  of  a  viscous  Rayleigh  instability  in  a  Stokes 
layer  profile  at  Re =200. 


Fig.  32.  Vorticity  contours  for  the  nonlinear  develop¬ 
ment  of  a  viscous  Rayleigh  instability  in  a  Stokes  layer 
profile  at  Re =200. 


71 


Fig.  33.  Wall  shear  stress  for  the  nonlinear  development 
of  a  Rayleigh  instability  in  a  Stokes  layer  profile  at 
Re =200,  showing  soliton  like  behavior  at  large  time. 


Fig.  34.  Streamfunction  contours  for  the  nonlinear  de¬ 
velopment  of  a  viscous  Rayleigh  instability  in  a  Stokes 
layer  at  Re =400.  The  vorticity  contours  for  this  solution 
show  the  spurious  oscillations  normally  encountered  in 
the  central  difference  scheme  for  the  high  Re  viscous 
problems  and  the  inviscid  problem. 


Fig.  35.  Vorticity  contours  for  the  nonlinear  develop¬ 
ment  of  a  viscous  Rayleigh  instability  in  a  Stokes  layer  at 
Re =400,  showing  spurious  oscillations. 
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Fig.  36.  Typical  streamfunction  contours  for  the  nonlin¬ 
ear  development  of  an  inviscid  Rayleigh  instability  in  a 
Stokes  layer,  using  central  difference  scheme. 
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Fig.  37.  Wall  slip  velocity  for  the  nonlinear  development 
of  an  inviscid  Rayleigh  instability  in  a  Stokes  layer  profile 
using  a  central  difference  scheme. 
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Fig.  39.  Streamwise  velocity  profiles  at  a  typical  stream- 
wise  station  in  a  viscous  Rayleigh  instability. 
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Fig.  40.  Typical  mesh  used  for  the  Navier- Stokes  com¬ 
putations  on  a  parabola  at  angle  of  attack.  Up  to  500 
streamwise  points  have  been  used  in  the  clustered  region 
near  the  leading  edge. 
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Fig.  41.  Comparison  of  steady  state  limit  of  the  unsteady 
Navier— Stokes  solutions  with  Davis  (1972)  computa¬ 
tions  of  flow  past  parabola  at  zero  angle  of  attack. 
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Fig.  42.  Streamlines  for  Navier-Stokes  solution  past  pa¬ 
rabola  with  an  impulsively  applied  no— slip  condition. 
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Fig.  43.  Vorticity  contours  for  Navier-Stokes  solution 
past  parabola  with  an  impulsively  applied  no -slip  condi¬ 
tion. 
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Fig.  46.  Verification  of  streamwise  grid -independence 
using  wall  shear  for  Navier— Stokes  solution  past  parab¬ 
ola  undergoing  smooth  pitch -up. 
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Fig.  47.  Verification  of  normal  grid— independence  using 
wall  shear  for  Navier— Stokes  solution  past  parabola  un¬ 
dergoing  smooth  pitch-up. 
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Fig.  48.  Primary  and  secondary  separation  lines  for  Navi- 
er— Stokes  solution  past  parabola  undergoing  smooth 
pitch— up.  Also  shown  are  comparisons  at  different  grid 
spacings. 
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Fig.  49.  Navier— Stokes  solutions  for  an  impulsive 
change  in  angle  of  attack  which  initially  remains  at¬ 
tached.  Shown  are  local  disturbance  developments  and 
the  increase  in  growth  rate  as  the  base  flow  approaches 
separation. 
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Fig.  50.  Blow— up  of  the  local  oscillatory  shear  stress  in 
figure  49,  showing  the  multiple  separations. 
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Fig.  51.  Vorticity  contours  in  the  boundary  layer  corre¬ 
sponding  to  figure  50.  Navier— Stokes  solutions  for  an 
impulsive  change  in  angle  of  attack  which  initially  re¬ 
mains  attached. 
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Fig.  52.  Effect  of  surface  perturbations  on  Navier- 
Stokes  solutions  for  impulsive  changes  in  angle  of  attack 
in  attached  flows.  Disturbances  increase  final  distur¬ 
bance  amplitude. 
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Fig.  53.  Vorticity  contours  in  the  boundary  layer  corre¬ 
sponding  to  figure  52.  Navier- Stokes  solutions  for  an 
impulsive  change  in  angle  of  attack  which  initially  re¬ 
mains  attached,  with  surface  perturbations. 
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Fig.  54.  Streamfunction  contours  in  the  boundary  layer 
corresponding  to  figure  52.  Navier— Stokes  solutions  for 
an  impulsive  change  in  angle  of  attack  which  initially  re¬ 
mains  attached  (with  surface  perturbations). 
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Fig.  55.  Comparison  of  steady  state  limit  of  the  unsteady 
Navier— Stokes  CUD-3  solutions  with  Davis  (1972) 
computations  of  flow  past  parabola  at  zero  angle  of  at¬ 
tack. 
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Fig,  56.  Comparison  of  Navier-Stokes  CUD— 3  skin 
friction  with  central  difference  skin  friction  at  early  times 
in  a  parabola  undergoing  smooth  pitch— up. 


